REPORT  DOCUMENTATION  PAGE 


AFRL-SR-BL-TR-Ol- 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instrui 
data  needed,  and  completing  and  reviewing  this  collection  of  Information,  Send  comments  regarding  this  burden  estimate  or  any  other  aspi 
this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-^1 88 

4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to 
valid  0MB  control  number,  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 

£>^Lf 

1 .  REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE 

-OT-  aool  Special  FINAL  Technical 

3.  DATES  COVERED  (From  -  To) 

15  Feb  00  -  30  Nov  00 

4.  TITLE  AND  SUBTITLE 

Conference  on  p  and  hp  finite  element  methods 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

F49620-00-1-0165 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

Manil  Suri,  P.I. 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

University  of  Maryland,  Baltimore  County 

1000  Hilltop  Circle 

Baltimore,  MD  21250 

8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 

9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Office  of  Scientific  Research 

801  North  Randolph  Street,  Room  732 

Arlington,  VA  22203-1977 

WRFORCEOFRCEC 

10.  SPONSOR/MONITOR’S  ACRONYM(S) 

AFOSR/NM 

11.  SPONSOR/MONITOR’S  REPORT 

Tfn  r\Tir  TtilfiTFCnKsCAU - 

12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT  tlOTlUtUf 

Distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

Enclosed  is  final  symposium  report  together  with  list  of  attendees, 
conditions . 

2001 

as  required  by  grant 

10817  050 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF: 


a.  REPORT 

Unclassified 


b.  ABSTRACT 

Unclassified 


c.  THIS  PAGE 

Unclassified 


17.  LIMITATION 
OF  ABSTRACT 


18.  NUMBER 
OF  PAGES 


^^3 


Manil  Suri 


19b.  TELEPHONE  NUMBER  (include  area 
code) 

(301)  588-3140 


Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  Z39.18 


p  and  hp  Finite  Element  Methods: 
Mathematics  and  Engineering  Practice 


Codfereiice  in  honor  of  the  65  th  birthday  of  Professor 
BarnaA.  Szabo 


Conference  Sponsored  in  Part  by: 

•  Air  Force  Office  of  Scientific  Research,  Air  Force 
Research  Laboratory,  USAF 

•  US  Army  Research  Office 

.  Mechanlcai  Eng.  and  System  Seienee  and  Math. 
Depts.  at  Washington  University. 


Summaries  of  Papers 

May  31  -  June  2, 2000 
Washington  University,  St.  Louis,  USA 


This  booklet  contains  the  abstracts  of  all  papers  presented  at  pFEM2000. 
The  abstracts  are  listed  in  alphabetical  order  of  first  author  of  each  paper, 
started  with  tiie  Plenary  Talks,  followed  by  the  Parallel  Talks.  For 
reader's  convenience,  the  abstract  number  is  enclosed  as  well. 

The  booklet  contans  three  parts: 

•  Table  of  Contents 

•  Completed  Abstracts 

•  Author  index  and  Addresses 


Selected  papers  presented  at  p-FEM2000  will  be  published  in  a  special 
issue  of  Ae  journal 

International  Journal  for  Numerical  Methods  in  Engineering 
and  a  special  issue  of  the  journal 

Computers  and  Mathematics  with  Applications 


Plenary  Lectures 


Borje  Andersson . 1 

Statistical  analysis  of  3D  multiple-site  fatigue  crack  growth 
in  aging  aircraft  using  a  hp-version  of  FEM 

I  VO  BabuSka . 3 

The  hp  version  of  fern.  From  birth  to  maturity  and  beyond 

Ted  Belytschko . 4 

Discontinuities  in  FEM  and  Meshfree  Methods  and  Level  Sets 

Benqi  Guo . 6 

Thirty  years  of  the  p  and  hp  versions  of  FEM  -  A  historical 
survey 

George  Karniadakis . . . 7 

Unstructured  spectral/hp  elements  for  turbulence  simulations 

Saeed  Paydarfar . . 8 

Why  2p  is  not  enough  for  Y2K  and  beyond 

Christoph  Schwab . 9 

Multiscale  high-order  FEM  for  plates  and  shells 

J.  Tinsley  Oden . . 10 

Use  of  hp  finite  element  methods  in  material  modeling 

Michael  J.  Wheeler . 11 

A  review  of  the  corporate  impact  of  P  technology  pioneered  by 
Prof.  Szabd 

Pax&llBl  Lectures 

Ricardo  Actis  and  Barna  Szab6 . 12 

Analysis  of  bonded  and  fastened  repairs  by  the  p-version  of  the 
finite  element  method  (AW  -  16) 

Slimane  Adjerid,  Joseph  Flaherty  and  Lilia  Krivodonova . 13 

A  posteriori  error  Estimation  for  Hyperbolic  Partial  Differential 
Equations  (A  -  51) 

Mark  Ainsworth . 14 

Some  remarks  on  a  posteriori  error  estimation  for  parameter 

dependent  probl&ns  (A  -  52) 

Mark  Ainsworth . 14 

The  stability  of  high  order  mixed  finite  element  families  for 
incompressible  flow  (A  -  53) 

Gal  Amar  and  Zohar  Yosibash . 15 

p-FEM  for  formulating  an  elastic  criterion  for  predicting 

mechanical  failures  at  2-D  singruJar  points  (GW  -  21) 

Franco  Auteri,  Luigi  Quartapelle  and  L.  Vigevano . 18 

Helmholtz— Neumann  spectral  solver  with  an  essential  treatment  of 
boundary  conditions  method  (A  -  56) 


i 


Franco  Auteri,  F,  Saleri  and  L.  Vigevano . 21 

Incompressible  Navier-Stokes  solutions  by  a  triangular  spectral/hp 
element  projection  method  (A  -  55) 

Franco  Auteri,  J.-L.  Guermond,  N.  Parolini  and  L.  Quartapelle . .24, 

On  the  relevance  of  the  inf -sup  condition  in  spectral  projection 
methods  (A  -  54) 


Mejdi  Azalez,  F.  Ben  Belgacem  and  M.  Karimi-Fard . 26 

Time  splitting/spectral  element  method  for  Navier-Stokes  equations 
(A  -  57) 

P.K.  Basu,  A.B.  Jorge,  S.  Badri  and  J.  Lin . ..28 

Higher  order  modeling  of  continue  by  finite  element^ 
boundaryelement,  and  wavelet  methods  (BW  -  1) 

Faker  Ben  Belgacem . 29 

The  hp — mortar  finite  element  method  for  mixed  elasticity  and 
stokes,  problems  (B  -  58) 

Edgar  Bertoti . . 


Dual-mixed  p  and  hp  finite  elements  for  elasticity  problems  (B  — 

59) 


Gabriella  Bognar  and  Tamas  Szab6 . 33 

The  solution  of  a  nonlinear  eigenvalue  problem  using  p-version  FEM 

(GW  >  2) 

Susanne  C.  Brenner . . 

Multigrid  methods  for  two-dimensional  interface  problems  (B  -  59a) 

Wei  Cai . . 

High  order  mixed  rwg  basis  functions  for  electromagnetic 
scattering  in  multilayered  media  and  applications  (C  -  60) 

T.  Cao  and  Don  W.  Kelly . 37 

Pointwise  and  local  error  estimates  for  the  quantities  of  interest 
in  two  dimensional  elasticity  (CW  -  24) 

Eduardo  Chan . . 

Interfacing  h-element  and  p-element  models  using  static 
condensation  (CW  -  3) 

Nikos  Chrisochoides  and  Demian  Nave . . 

Parallel  guaranteed-quality  h-refinement  and  mesh  generation  (C  - 

60a) 

Fr4d4ric  Cugnon . . 

Error  estimation,  automatic  p-  and  hp-adaptation  using 
hierarchical  concept  (CW-25) 

Leszek  Demkowicz . . . . . . 43 

Adaptive  hp-FE  modeling  for  maxwell 's  equations  (C  -  61) 


f 


t 


ii 


Karen  Devine,  Bruce  Hendrickson,  Erik  Boman,  Courtenay  Vaughan  , and 

Matthew  St.  John . . . 45 

Zoltan:  A  library  for  dynamic  load  balancing  (D  -  62) 

Saikat  Dey,  J.  J.  Shirron  and  L.  S.  Couchman . 47 

Parameterized  modeling  of  multi-layered  stiffened  shells  for  p- ' 
version  analysis  (D  -  63) 

Stephan  P.  Engelstad . 48 

Development  of  p-version  handbook  solutions  for  analysis  of 
composite  bonded  joints  (EW  -  21) 

Joseph  E.  Flaherty  and  Jean-Francois  Remade . 50 

Adaptiave  discontinuous  Galerkin  method  with  orthonormal  basis  (F 

-  63a) 

Joseph  E.  Flaherty,  Mark  S.  Shephard  and  J.  D.  Teresco . 51 

Architecture-dependent  finite  element  computation  using  the 
Rensselaer  partition  model  (F  -  64) 

J.  de  Frutos  and  Julia  Novo . 54 

A  posteriori  error  estimation  for  evolutionary  disipative 
equations  (F  -  65) 

Tat  Ching  Fung  and  Sieng  Khuan  Chow . 56 

Higher  order  accurate  complex  time  step  mmethod  (FW  -  4) 

Klaus  Gerdes,  Malkus  Melenk  and  Christoph  Schwab . . . 57 

Fully  discrete  hp  finite  elements  (6  -  66) 

Ben-yu  Guo . 58 

Jacobi  spectral  method  <6  -  68) 

Xian-Zhong  Guo,  I.  Norman  Katz  and  Ning  Hu . 61 

Multi-p  processes  and  pre-conditioners  (XW  -  19) 

Marc  Halpern . 62 

P-adaptive  technology:  perspectives  on  its  impact  and  future 
growth  (HW  -  5) 

Herbert  Heuer . 64 

Preconditioned  iterative  solvers  for  coupled  FEM-BEM  systems  (H  - 

68a) 

Stephan  Holzer . . : . 65 

On  nonlinear  and  higher-order  plate  analysis  using  the  p-version 
of  the  FEM  (HW  -  6) 

Endel  V.  larve . 66 

Spline  variational  three  dimensional  stress  analysis  of  laminated 
composite  plates  with  open  holes  (IW  -  23) 

Tim  W.  Jackson . 67 

Practical  application  of  high-order  methods  in  engineering 
practice  (JW  -  7) 

David  A.  Kopriva,  Stephen  L.  Woodruff  and  M.Y.  Hussaini . 68 

Computation  of  electromagnetic  scattering  with  a  non- 
conformingdiscontinuous  spectral  element  method  (K  -  69) 


iii 


Roland  Krause  and  Dominik  Scholz . 70, 

On  the  application  of  high  order  methods  to  problems  in  structural 
dynamics  (K  -  70) 

Ric  Leist  and  Yupu  Shi . 71 

The  application  of  p-Version  finite  element  solutions  for  advanced  - 
fatigue  analysis  (LW  •-  8) 

Matthias  Maischak . 72 

Multiplicative  schwarz  algorithms  for  the  Galerkin  boundary 
element  method  of  the  first  kind  (M  -  70a) 

Edward  A.  W.  Maunder . 73 

Non  conventional  dual  analyses  in  the  context  of  error  estimation 

(MM  -  9) 

Rainer  Niekanqp  and  Erwin  Stein . 76 

Hierarchical  multilevel  graded  hp-FE  algorithms  with 
parallelprocessing  based  on  hexaeder  elements  with  node  regular 
refinements  (M  -  71) 

Hae-Soo  Oh . . 

Highly  accurate  mode-separated  energy  release  rates  for 
delaminated  composite  laminates  (O  -  72) 

Stephan  Ohnimus . . 

Equilibrated  residual  spatial  h-  and  p-  error  estimators  for 
elliptic,  parabolic  and  elastoplastic  problems  (O  -  73) 

Istvan  Paczelt  and  Tamas  Szab6 . 81 

Solution  of  the  contact  optimization  problems  of  cylindrical 
bodies  using  hp-FEM  (PW  -  10) 

Abani  Patra,  J.  Long,  A.  Bauer . 83 

On  the  development  of  large  scale  parallel  adaptive  3D  codes:  data 
management  schemes  and  multi-level  iterative  substructuring 
solvers  (P  -  73a) 

Scott  Prost-Domasky,  C.  Brooks  and  K.  Honeycutt . 85 

The  application  of  p-version  finite  element  methods  to  fracture- 
dominated  problems  encountered  in  engineering  practice  (PW  -  11) 

Sivapriya  Ramanathan  and  Manish  Parashar . 86 

Adaptive  distribution  of  dynamic  grid  hierarchies  (R  -  73b) 


Ernst  Rank,  Henrike  Broker,  and  Alexander  DOster . 88 

The  p-version  for  linear  and  non-linear  analysis  in  three 
dimensions  (BST  -  14) 

Paul  Reid . . 89 

Current  usage  of  two  p-version  finite  element  analysis  codes  and 
proposed  concepts  for  future  growth  and  capability  (RM  -  15) 

Beatrice  Riviere  and  Mary  F.  Wheeler . 90 


A  posteriori  error  estimates  for  a  discontinuous  galerkin  method 
applied  to  elliptic  problems  (R  -  74) 


f 


f 


IV 


Anna-Margarete  Saendig . .  91 

Fem  -  error  estimates  for  stationary  Navier-Stoke  equations  in 
nonsmooth  domains  (S  -  75) 


John  Schiermeier,  R.  K.  Kansakar,  D.  Mong  ,  J.  B  .Ransom  ,  M.  A. 

Aminpour  ,  W.  J.  Stroud . 92 

p-Version  interface  elements  in  global/local  analysis  (SW  -  17) 

Kirk  Schloegel,  George  Karypis  and  Vipin  Kumar . 94 

Parallel  and  adaptive  multi-constraint  graph  partitioning  <S  -75a) 

Oominik  Schoetzau  and  Christoph  Schwab . 96 

The  hp-version  of  the  discontinuous  galerkin  time-stepping  method 

(S  -  76) 

H.  Schulz  and  Wolfgang  L.  Wendland . 97 

On  a  high  order  boundary  element  method  (S  -  76a) 

Padmanabhan  Seshaiyer . 100 

Non-Conforming  hp  mortar  finite  .element  methods..  {3  -  77) 

Jie  Shen.  . . . 101 

Stable  and  efficient  spectral  methods  in  unbounded  domains  using 
laguerre  functions  (8  -  78) 

Mark  S.  Shephard,  R.  M.  O'Bara,  X.  Luo  and  X.  Li . 103 

Algorithm  to  generate  p-version  meshes  for  curved  domains  (SW  - 

18) 


Mark  S.  Shephard,  Joseph  E.  Flaherty,  Ottmar  Klaas,  Jim  D.  Teresco 

and  Wes  Turner . 105 

Trellis  -  A  geometry  based  adaptive  analysis  framework  for  a 
multi-processor  environment  (SW  -  20) 

Spencer  J.  Sherwin  and  Joaquim  Peiro . 107 

Mesh  adaption  for  improved  geometry  representation  with  high  order 
elements  (S  -  79) 

Spencer  J.  Sherwin,  Peiro  Joaquim  &  Doorly  D.J . 109 

Computational  haemodynamics  in  arterial  bypass  grafts  (S  -  80) 

Ernst  P.  Stephan... . Ill 

Additive  Schwarz  methods  for  the  hp-version  of  the  boundary 
element  method  for  first  kind  integral  equations  in  R3  (S  -  80a) 

Manil  Suri . 112 

Reliability  of  an  hp  algorithm  for  buckling  analysis  (S  -  81) 


Yehuda  Volpert  and  Scott  Prost-Domasky . 114 

The  challenge  of  obtaining  reliable  results  in  contact  problems 

(VW  “  12) 

Tim  Warburton  and  Jan  S.  Hesthaven . 116 

Application  of  nodal  unstructured  spectral  elements  to  Maxwell's 
equations  (W  -  82) 


f 


f 


V 


Beth  Wingate  and  Mark  Taylor . 

Diagonal  mass  matrix  spectral 

(W  -  83) 


. . .118 

methods  for  triangles  and  simplices 


Zohar  Yosibash . 121 

Extracting  edge  flux  intensity  functions  for  the  laplacian  using 
p-FEM  (Y  -  83a) 

Lily  Zhang  and  Barna  Szab6 . 122 

The  solution  of  viscoelastic  probl&as  by  the  mixed  hp-version  of 
the  finite  element  method  (Z  84) 


vi 


STATISTICAL  ANALYSIS  OF  3D  MULTIPLE-SITE 
FATIGUE  CRACK  GROWTH  IN  AGING  AIRCRAFT 
USING  A  VERSION  OF  FEM 


Borje  Andersson 

The  Aeronautical  Research  Institute  of  Sweden  Box  11021,  S-161 11  Bromma, 

Sweden  (ba@ffa.se) 


The  present  paper  deals  with  statistical  analysis  of  fatigue  crack  growth  and  crack 
coalescence  in  riveted  joints  in  real-life  aircraft  structures.  The  practical  problem 
in  minH  is  a  part  of  the  “Ageing  Aircraft”  problem,  i.e.  the  case  when  old  aircraft 
exhibit  many  small  fatigue  cracks  which,  in  relatively  few  load  cycles,  link-up  to 
form  large  cracks  which  jeopardize  aircraft  safety. 

In  the  mathematical  model,  several  simplifications  are  introduced  what  regards  con¬ 
tact  and  friction  between  rivets/skin  and  skin/stiffeners.  The  equations  of  3D  elas¬ 
ticity  ^6  assumed  valid,  i.e.  in  standard  notation, 

Cijki  Uj,ki{t)  =0  on  Q 

(1)  Cijki  nj  Uk/t)  =  Ti^(l  +  sin(cjt))  on  ^r(t,K) 

«i(t)  =0  on 

where  u  =  (tii,U2,«3)^  is  the  displacement  vector,  T  =  (Ti,r2,T3)^  the  traction 
vector,  n  =  (ni,n2,n3)^  the  normal  vector,  t  is  time,  w  an  angular  frequency  such 
that  inertia  effects  become  negligible.  The  domsun  with  surface  r(t)  =  ■*'^r(t)U®r, 
ffnnt.a.ins  K{t)  C  I  cracks  of  sizes  and  shapes  characterised  by  a  few  parameters 
{oi(t),5<(t),r<(t)  1  1  <  t  <  iif(t)}.  Crack  sizes  increase  with  time  due  to  the  time 
varying  load  (fatigue).  The  domain  of  interest  is  that  of  a  typical  aircraft  fuselage 
with  spars,  frames,  doublers  and  rivets,  all  parts  modelled  as  full  3D  objects 
in  the  numerical  analysis.  This  leads,  in  the  F&analysis  to  algebraic  systems  of 
equations  having  >  10®  degrees  of  freedom.  In  an  experimental  study  it  was  found 
that  the  model  used  accurately  predicts  multiple-site  fatigue  crack  growth  in  a  butt 
joint,  given  the  initial  crack  pattern. 

A  Monte  Carlo  giTniilfl.tinn  technique  is  used  to  calculate  statistical  distributions  of 
the  fatigue  life  for  thousands  of  different  initial  crack  flaw  distributions  in  aircraft. 
The  fatigue  crack  growth,  crack  coalescence,  crack  branching  (“flapping”)  is  fol¬ 
lowed  deterministically,  with  control  of  error,  for  each  of  the  initial  damage  scenario 
considered. 

This  type  of  analysis  requires  typically  10®-10^  accurate  solutions  of  (1)  for  various 
boundaries  ^r(t,  if)  (crack  sizes/distributions).  A  mathematical  splitting  method, 
which  provide  the  high  efficiency  needed,  is  used.  In  the  splitting  scheme,  the 
problem  is  divided  into  subproblems,  which  are  solved  using  a  hp- version  of  the  finite 
element  method,  where  the  solution  sought  is  obtained  by  proper  superposition. 


Convergence  to  exact  solution  for  stress  intensity  factors  is  exponential  in  the  pre- 
asymptotic  range. 

For  statistical  analysis  on  domains  of  present  complexity,  the  computational  effort 
might  be  very  large.  However,  the  mathematical  splitting  of  the  problem  can  be 
mapped  very  efficiently  on  a  cluster  of  SMP-computers.  A  new  solution  paradigm 
is  developed  for  control  of  the  analysis  which  shows  excellent  scalahUity.  On  a  60 
server  cluster  having  in  all  1000  CPu’s,  a  computational  speed  >  0.5  TFLOPS  was 
recorded.  Finally,  results  from  large-scale  analysis  of  a  large  stiffened  panel  are 
presented  and  conclusions  drawn. 


THE  p- VERSION  OF  FEM,  FROM  BIRTH  TO 
MATURITY  AND  BEYOND 


Ivo  BabuSka 

Texas  Tnstit.iit.ft  for  Computational  and  Applied  Mathematics  (TIC AM),  The 
University  of  Texas  at  Austin,  SHC  304, 105  W.  26th  St.,  Austin,  TX  78712,  USA 

(buska@ticam.utexas.edu) . 


The  first  part  of  the  talk  will  briefiy  address  the  ideas  behind  the  first  phase  of  the 
development  of  the  p- Version,  especially  with  respect  to  the  work  of  B.  Szabo  and 
the  major  results  obtained  later.  The  second,  part  of.  the  lecture  will  elaborate  on 
the  generalized  finite  element  method  (GFEM)  which  uses  special  shape  functions 
related  to  the  problem.  This  leads  to  the  higher  accuracy  of  the  solution  and  its 
efficiency.  GFEM  can  also  be  used  to  avoid  the  need  of  the  mesh  generator  for  very 
complex  domains  for  example  when  the  domain  has  many  inclusions.  GFEM  is  very 
fiexible  and  utilizes  all  advantages  of  the  standard  h-  and  p-  Version  of  FEM. 
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DISCONTINUTTES  IN  FINITE  ELEMENTS  AND 
MESHFREE  METHODS  AND  LEVEL  SETS 

Ted  Betytschko 


Mechanical  EngiiKering,  Northwestern  University,  Evanston,  IL  60208  USA 

(tedbelytschko@nwu.edu) 


This  p^r  presents  new  methods  for  the  construction  of  functions  with  arbitrary 
discontinuities  and  discontinuous  derivatives.  In  finite  element  methods,  these 
discontinuities  are  completely  independent  of  the  finite  element  mesh.  In  meshless 
methods,  they  have  advantages  over  the  visibility  techniques  developed  in  Belytschko,  Lu 
and  Gu  [1].  The  techniques  are  based  on  the  concepts  originated  in  Belytschko  and  Black 
[2]  and  Moen  et  al.  [3].  The  applications  of  the  techniques  are  diverse:  cracks  and  crack 
growth,  shw  bands,  and  phase  changes  are  some  of  the  areas  of  modeling  which  require 
discontinuities.  Although  evolving  discontinuities  are  often  treated  by  remeshing,  this 
becomes  quite  awkward  when  the  evolution  of  many  discontinuities  are  considered. 

Let  the  discontinuities  in  tire  dependent  variable  u(x)  occur  on  surfia^esr^,  vdtich 
are  described  by  signed  distance  functions  /(x)=  minflx-'^j/gn(n-(x-I)),  where  i  is 

a  point  on  the  surface  of  discontinuity  r„  and  n  is  a  unit  normal  to  the  surfiice  of 
discontinuity.  The  discontinuity  surface  is  represented  in  the  meshfree  and  finite  element 
models  by  /(x)  =^JV/(x),  vdiere  repeated  indices  herein  are  summed  over  the 
appropriate  range.  The  position  of  the  discontinuity  can  be  updated  by  level  set  methods. 
When  the  su[^ort  is  bisected,  the  iqrproximation  is  given  by 

“(X) = Z  + a/ffO;  W)) 

where /r(x)  is  a  modified  step  function.  The  coefficient  of  aj  is  called  the  enrichment 
For  a  slit  support,  i.e.  a  support  in  which  the  discontinuity  ends,  the  qrproximation,  is 


where  bfl(x)  are  branch  function  around  the  discontinuity.  Generally  more  than  one 
branch  function  is  needed  ftir  each  enriched  node;  examples  are 

^i(*)“|y‘sm^,r^sin— J  or  6jf(x)=|V^sin— ,V^sinjsin^,Vrcosj,Vrcosjsin^J 


As  can  be  seeri  in  the  branch  functions  on  the  left,  th^  are  discontinuous  across  the  line 
/r  but  continuous  and  well  behaved  in  the  domain  surrounding  the  discontinuity.  For 
linear  fracture  mechanics,  the  branch  functions  on  the  right  have  been  used  in  Belytschko 
and  Black  [2].  The  last  three  functions  are  not  discontinuous,  but  were  added  to  speed  the 
corivergence  to  elastic  fincture  problems.  Enrichments  for  discontinuities  in  the 
derivatives  and  for  single  components  of  a  vector  function  can  be  constructed  similarly. 
The  latter  can  be  used  to  model  shear  bands.  Enrichments  have  also  been  developed  for 
intersecting  discontinuities,  see  [4]. 
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THIRTY  YEARS  OF  THE  p  AND  h-p  FINITE  ELEMENT 
METHODS  -  A  HISTORICAL  SURVEY 

Benqi  Guo 

Department  of  Mathematics,  University  of  Manitoba,  Canada 
(bguo®  newton.amath.  umanitoba) 


As  two  major  approaches  of  today’s  finite  element  method,  the  p  and  h-p  versions 
have  significantly  developed  in  the  past  three  decades.  The  development  of  the  p  and 
h-p  finite  element  methods  not  only  provides  a  reliable  computational  tool  to  modem 
engineering/scientific  computing,,  but  also  offers  fresh  insights  into  many  engineering 
and  scientific  problems.  It  has  substantially  enriched  the  knowledge  of  the  finite 
element  method  and  many  relevant  fields  of  mathematics  such  as  approximation 
theory,  regularity  theory  of  partial  differential  equations  with  piecewise  analytical 
data,  and  theory  of  shells  and  plates. 

Originated  from  the  computation  of  stroctural  mechanical  problems,  the  p  and  h-p 
finite  element  methods  have  made  a  splendid  journey,  from  the  birth  in  the  later 
1960’s  to  today’s  state  of  the  art.  The  p  and  h-p  finite  element  methods  are  now 
widely  used  in  almost  every  branch  of  engineering.  Various  large-scale  minnriprriai 
and  research  codes  of  the  p  and  h-p  version  are  available  to  engineering  computa¬ 
tion  in  industry  as  well  as  many  research  institutes  and  universities.  These  codes 
furnished  with  well-developed  algorithms  and  combined  with  modem  computer  tech¬ 
nology  make  it  possible  to  solve  engineering  problems  of  great  complexity.  The  p 
and  h-p  finite  element  methods  have  successfully  penetrated  into  today’s  engineering 
practices  and  academic  research. 

The  research  on  the  theory  of  the  p  and  h-p  finite  element  method,  started  in  the 
later  1970  s,  has  established  a  strong  mathematical  foundation,  which  is  able  to 
sufficiently  support  the  scientific  and  engineering  computing  and  to  offer  new  ways 
of  improving  the  effectiveness  of  engineering  computations. 

The  tremendous  efforts  and  contributions  from  engineers  and  mathematicians  lead 
to  the  great  achievements  of  the  p  and  h-p  finite  element  methods  in  the  twentieth 
century.  In  the  year  of  2000,  it  is  worth  while  to  address  the  history  of  the  p  and  h-p 
finite  element  methods  for  better  understanding  in  the  future.  The  survey  will  give 
a  review  of  major  developments  of  the  p  and  h-p  finite  element  methods  in  past  30 
years  on  the  following  aspects:  approximation,  algorithms,  engineering  applications, 
mathematical  framework,  and  solvers  of  linear  algebraic  equations. 

The  survey  will  address  some  outstanding  problems  of  p  and  h-p  finite  element 
methods  in  the  beginning  of  new  millftnniiiTn 
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UNSTRUCTURED  SPECTRAL/hp  ELEMENTS  FOR 
TURBULENCE  SIMULATIONS 


George  Karniadakis 

Division  of  Applied  Mathematics,  Brown  University,  Providence,  RI  02912,  USA 

(gk@c£m.brown.edu) 


In  the  presentation  we  will  discuss  several  versions  of  high-order  bases  (some  inspired 
by  Prof.  Szabo’s  pioneering  work)  and  present  algorithms  and  results  for  simula¬ 
tion  of  turbulent  flows  in  complex-geometry  domains.  In  particular,  new  results  on 
turbulent  wakes  at  Reynolds  number  of  100,000  will  be  presented  for  first  time. 
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WHY  2p  IS  NOT  ENOUGH  FOR  Y2K  AND  BEYOND 

Saeed  Paydarfar 


Vice  President  and  General  Manager,  Advanced  Enterprise  Solutions,  CA,  USA 
(spaydarfar@aes4solutions.com) 


Today,  we  arc  ^ly  in  the  midst  of  another  Industrial  Revolution:  a  digital  revolution, 
i^d,  just  as  m  the  case  of  die  first,  those  oiganizations  that  can  quickly  attain  mastery  of 
the  new  ^igm  will  excel  beyond  their  competition,  and  those  that  don’t  will  quickly 
^  th«  there  ^  schedule  and  quality  will  not  meet  the  new  standards  of  the  time. 

of  the  first,  rapid  deployment  of  new  technology  is  met  with 
sigrafic^  techmc^  and  cultural  challenges,  often  defeating  most  well  intentioned  efforts 
to  provide  newels  and  integration  solutions.  This-industriaLrevolution  spans  the  entire 
^hn^ultural  foundation  of  our  everyday  life  effecting  everything  fiom  ^  we  buy  a 
book  to  how  we  develop  advanced  products  in  a  fully  expanded  digital  enterprise. 

^e  yw  2000  ^2K)  and  beyond  will  mark  toe  beginning  of  a  new  phase  of  the  digital 

into  the°hi  h”  decades,  software  developers  have  put  enormous  capability 

into  toe  hw^  of  the  industrial  community  for  the  development  of  pnxlucts.  T^y,  tte 

challenge  is  how  to  mtegrate  these  software  tools  and,  subsequentiy,  what  appears  to  be 

cmfiguration  comnillcd  envitonniem  that  allows  the  seemless  flow  of 
«^on  ^  ^  discipitoe  to  the  next  New  o«h„ds  thHow 

assurance,  and  increased  functional  performance  are  the 


!?™te  Element  ai^ysis  is  a  key  component  of  toe  product  development  process 
^velop^nte  in  FE  methods  have  certainly  provided  enormous  benefit^to  toe  Analysis 
wnmunity.  In  rer^t  times,  of  particular  value  is  the  higher  order  polynomial,  p-versfon, 
method.  This  method  has  a  robust  integration  with  toe  design  geometry 
^  Its  bound^  conditions,  provides  an  error  estimate  of  the  analysis,  and  proviS 
functional  methods  that  are  not  easily  achieved  in  the  tradition  FE  ^proach. 


^wi?  to  enterprise  wide  integration  and  toe  enabling 

^  requiTO  a  stall  mrx  of  enormous  proportion.  Tool  benchmarking  of  toe  pal 
ne^  to  be  r^laced  by  mtegration  projects  that  focus  on  culture  and  global  intraration. 
Integration  roluhons  need  to  be  enterprise  wide,  provide  stable  functionality,  crafted  for 

contimmus  improvement,  and  provide  methods  for  an  older  of  magnitude  reduction  in 
cost  and  schedule  with  significantly  increased  quahty.  reuucnon  rn 
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MULTISCALE  HIGH  ORDER  FEM  FOR  PLATES  AND 

SHELLS 


Christoph  Schwab 

Seminar  fur  Angewandte  Mathematik,  ETH  Zurich,  Ramistrasse  101,  CH8092 
Zurich,  Switzerland  (schwab@sam.math.ethz.ch) 


The  efficient  numerical  modelling  of  laminated  composites  has  received  increasing 
attention  in  recent  years  -  examples  are  sandwich  plates  and  shells,  Hber-reinforced 
composites  and  the  like. 

While  global  response  to  external  loadings  can  be  reliably  assessed  on  the  basis  of 
averag^,  or  homogenized,  models,  interlamina  stresses  can  only  be  obtained  by 
resolution  of  the  small  scales  of  the  problem. 

This  requirement  of  SCALE  RESOLUTION  and  the  low  solution  regularity  at  the 
matrix-composite  interface  seems  to  conffict  with  the  use  of  high  order,  spectral  or 
hp-FE  approaches  which  is  commonly  based  on  large  elements  with  smooth  shape 
functions  of  high  order. 

We  present  a  new  approach  for  the  design  of  high  order  multiscale  FEM  for  sandwich 
and  lattice  materials.  The  elements  discretize  concurrently  macro-  and  microstruc¬ 
ture  of  the  material. 

The  FEM  realize  exponential  convergence  independently  of  the  number  of  layers  in 
the  sandwich/  the  number  of  fibers  per  unit  length.  It  is  based  on  a  new  spectral 
approach  to  homogenization  and  on  nonpolynomial  shape  (or  director)  functions  in 
the  Finite  Element  Analysis. 

Numerical  results  for  composite  material  with  periodic  structure  and  for  sandwich 
plates  and  shells  will  be  presented.  This  is  joint  work  with  A.M.  Matache  (student), 
I.  Babuska  (Texas),  B.A.  Szabo  and  R.  Actis  (St.  Louis,  USA). 


USE  OF  hp  FINITE  ELEMENT  METHODS  IN 
MATERIAL  MODELING 


J.  Tinsley  Oden 

The  Texas  Institute  for  Computational  and  Applied  Mathematics,  The  University 
of  Texas  at  Austin,  105  W.  25th  Street,  SHC  304,  Austin,  Texas  78712,  USA 

(oden@ticam.utexas.edu) 


In  this  lecture,  we  discuss  several  dif5culties  in  the  implementation  of  hp  and  p 
finite  element  methods.  We  describe  a  number  of  techniques  for  overcoming  these 
difficulties.  The  goal  is  to  develop  a  convenient  data  structure,  useful  for  developing 
robust  application  codes,  which  allows  the  implementation  of  mesh  adaptivity  with 
non-uniform  distributions  of  polynomial  order  of  p.  These  processes  involve  parallel 
nodal  constraint  algorithms  introduced  in  the  80s  [1],  the  use  of  PUMs  as  a  basis  of 
hp  approximations,  and  recent  implementations  of  discontinuous  Galerkin  methods. 

Next,  we  focus  on  applications  of  hp  methods  to  the  analysis  of  highly  heterogeneous 
materials.  We  demonstrate  how  parallel  hp  methods  can  be  used  as  an  analysis  tool 
for  adaptive  modeling,  as  well  as  adaptive  meshing.  Error  Estimations  and  are  also 
presented. 

Reference 

[1]  Oden,  J.T.,  Demkowicz,  L.F.,  Rachowicz,  W.,  and  Hardy,  O.  ’’Toward  a  Universal 
h-p  Adaptive  Finite  Element  Strategy”,  Computer  Methods  in  Applied  Mechanics 
and  Engineering,  three  parts:  77  (1989),  pp.  79-112,  pp.  113-180,  pp.  181-212. 
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A  REVIEW  OF  THE  CORPORATE  INPACT  OF 
p-TECHNOLOGY  PIONEERED  BY  PROF.  S2:AB0 

MichaelJ.  Wheeler 

Vice  President-Web  Stiate^es,  TechNet  International  Corp.,  Pittsburgh  PA,  USA 

(m\y*eeler@nauticom.net) 

In  the  early  1970’s,  Professor  Bama  Szabo  pioneered  an  exciting  new  way  to  enhaiice 
the  reliability  of  numerical  simulation.  This  has  given  rise  to  a  number  of  cotnpanies 
that  were  founded  to  commercialize  and  expand  this  technology  in  a  number  of  different 
areas.  This  presentation  will  review  these  corporations  and  their  contributions  to  the 
overall  technology  field. 
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ANALYSIS  OF  BONDED  AND  FASTENED  REPAIRS  BY 
THE  P-VERSION  OF  THE  FINITE  ELEMENT  METHOD 

Ricardo  L.  Actis 

Engineering  Software  Research  and  Development,  10845  Olive  Blvd.,  Suite  170,  St 
Louis,  MO  63141,  USA  (ricardo@esrd.com) 


Bama  A.  Szabo 

Center  for  Computational  Mechanics,  Washington  University,  Campus  Box  1129  St 
Louis,  MO  63130-4899,  USA  (szabo@ccm.wustl.edu) 


The  development  of  efficient  and  reliable  methods  for  the  analysis  of  fastened  structural 
connections  and  bonded  joints  is  among  the  most  important  problems  in  aerospace 
applications  because  these  connections  are  common  sites  of  failure  initiation. 


Re^istic  representations  of  festened  connections  must  account  for  the  diameter  of 
fetener  holes,  the  stiffiiess  of  the  fasteners,  the  effects  of  interference  fitting  or  looseness 
of  fit,  contact  ^tween  the  fasteners  and  the  connected  plates,  and  yielding  of  the  plates 
Co^uently  the  complexity  of  the  models  is  such  that  it  is  difficult  and  time  consuming 
to  treat  them  using  conventional  finite  element  or  boundary  element  procedures 


Similarly,  the  treatment  of  bonded  joints  is  complicated  by  the  fact  that  bonded  patehes 
are  t)pically  compos^  of  several  thin  layers  of  orthotropic  material  and  the  adhesive 
nidtcn&l  mEy  3^6ld  under  lofld,  recjuiring  e  nonlinesr  EnElysis. 


Fortunately,  these  problems  lend  themselves  to  standardization  in  FEA-based  handbook 
Iibranes.  Struc^  connections  can  be  described  in  a  parametric  form.  Therefore  a 
topolopcally  similar  femily  of  connections  need  to  be  meshed  only  once.  The  FEA- 
bas^  handbook  library  provides  for  the  archival  and  recall  of  entries  for  analysis  Users 
n^  only  enter  the  parameters.  The  finite  element  mesh,  being  associative,  is  adjusted 

autom^cally.  Therefore  users  need  not  be  concerned  with  meshing  or  other  details  of  the 
analysis  process. 


The  firate  element  solution  and  the  computation  of  all  data  of  interest  are  performed 
automatically  and  the  results  are  displayed  in  tabular  and  graphical  forms.  Automatic 
quality  control  tests  are  performed  and  reported  for  purposes  of  documentation.  There  are 
provisions  for  the  computation  of  margins  of  safety.  The  advanced  FEA-based  handbook 
librEiy  produces  results  of  sufficient  reliEbilily  to  permit  certificEtion  of  repEirs  End  joints 
on  the  basis  of  computed  information.  This  will  be  demonstrated  by  examples. 
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A  POSTERIORI  ERROR  ESTIMATION  FOR 
HYPERBOLIC  PARTIAL  DIFFERENTIAL  EQUATIONS 


Slimane  Adjerid 

Department  of  Mathematics,  Virginia  Polytechnic  Institute,  Bladcsburg,  VA 
24061,  USA  (adjerids@math.vt.edu) 


Joseph  E.  Flaherty  and  Lilia  Krivodonova 

Department  of  Computer  Science,  Rensselaer-,  Polytechmc  lnstitute,  Troy,  NY 
12180-3590,  USA  (flaherje/krilole@cs.rpi.edu) 


The  discontinuous  Galerkin  method  (DGM)  is  an  appeding  appro^  to  addr^ 
problems  having  discontinuities,  sudi  as  those  that  arise  in  hyperbolic  conservation 
laws.  Originally  developed  for  neutron  transport  problems,  the  DGM  has  been 
used  to  solve  ordinary  differential  equations  and  hyperbolic  ,  parabolic,  and  elliptic 
partial  differential  equations.  The  DGM  may  be  regarded  as  a  way  of  extending 
finite  volume  methods  to  arbitrarily  high  orders  of  accuracy.  The  basis  is  piecewise 
continuous  relative  to  a  structured  or  unstructured  mesh.  As  such,  it  can  sharply 
capture  discontinuities  in  the  solution.  Regardless  of  order,  the  DGM  has  a  simple 
communication  pattern  to  elements  with  a  common  face  that  makes  it  useful  for 
parallel  computation.  It  can  handle  problems  in  complex  geometries  to  high  order. 
In  order  for  the  DGM  to  be  useful  in  an  adaptive  setting  with  ^refinement  (mesh 
refinement  and  coarsening)  or  p-refinement  (method  order  variation),  techniques 
for  estimating  the  discretization  errors  should  be  available  both  to  gmde  adaptive 
enrichment  and  to  provide  a  stopping  criteria  for  the  solution  process.  We  wiU  show 
that  the  DGM  discretization  finite  element  solution  of  one-dimensional  hyperbolic 
conservation  laws  with  degree  p  exhibit  superconvergence  at  Radau  points  of  de^ee 
p  -H  1  with  the  fixed  endpoint  selected  at  the  upwind  end  of  each  element.  We 
use  this  superconvergence  result  to  construct  asymptotically  exact  a  posteriori  error 
estimates  for  hyperboHc  systems  of  partial  differential  equations.  Finally,  we  prraent 
numerical  results  for  several  computational  examples  that  show  the  effiaency  of  our 
a  posteriori  error  estimator. 
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SOME  REMARKS  ON  A  POSTERIORI  ERROR 
ESTIMATION  FOR  PARAMETER  DEPENDENT 

PROBLEMS 


Mark  Ainsworth 

Mathematics  Department,  Strathclyde  University,  Livingstone  Tower,  26 
Richmond  Street,  Glasgow  G1  IXH,  Scotland  (M.Ainsworth@strath.ac.uk) 


Standard  a  posteriori  error  estimators  when  applied  to-problems  containing  a  pa¬ 
rameter  are  often  found  to  be  sensitive  to  the  parameter  in  various  limitiTig  cases. 
This  either  means  that  as  the  parameter  approaches  certain  limits,  the  estimator 
will  either  become  increasingly  pessimistic  or  increasingly  unreliable.  The  talk  will 
present  examples  of  this  type  of  behaviour  and  discuss  methods  whereby  the  robust¬ 
ness  and  reliability  of  the  estimators  may  be  restored. 
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THE  STABILITY  OF  HIGH  ORDER  MIXED  FINITE 
ELEMENT  FAMILIES  FOR  INCOMPRESSIBLE  FLOW 

Mark  Ainsworth 

Mathematics  Department,  Strathclyde  University,  Livingstone  Tower,  26 
Richmond  Street,  Glasgow  G1  IXH,  Scotland  (M.Ainsworth@strath.ac.uk) 


The  talk  will  discuss  the  stability  of  mixed  finite  element  feuiiilies  on  meshes  contain¬ 
ing  high  aspect  ratio  aflBne  quadrilateral  elements.  Sharp  theoretical  estimates  for 
the  dependence  of  the  inf-sup  constant  on  the  aspect  ratio,  the  element  size  and  the 
spectral  order  of  the  elements  will  be  presented,  along  with  supporting  numerical 
evidence. 
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p-FEM  FOR  FORMULATING  AN  ELASTIC  CRITERION 
FOR  PREDICTING  MECHANICAL  FAILURES  AT  2  -D 

SINGULAR  POINTS 

Gal  Amar  and  Zohar  Yosibash 

Mechanical  Engiiieering  DepL  Ben-Gurion  University  of  the  Negev,  Beer-Sheva,  Israel 
(galam@bguinail.bgu.ac.il,  zohaiy@pversion.bgu.ac.il) 

Failures  of  mechanical  components  and  electronic  devices  usually  start  at 
reentrant  comers  or  the  intersection  of  multi-material  inter&:es.  To  date  the  structural 
engneering  community  is  lacking  a  criterion  for  predicting  the  failure  onset  in  the 
vicinity  of  these  singular  points.  In  linear  elastic  fiacture  mechanics  for  brittle  niaf«.riaic  |t 
is  well  known  that  a  single  parameter  of  the  solution  in  the  vicinity  of  the  crack  tip 
(called  stress  intensity  factor- JIT;)  can  be  correlated  to  a  material  depeiKlent  value  for 
predi^ng  fracture  onset.  Same  methodology  may  be  t^rplicable  for  singularities 
associated  with  sharp  re-entrant  comers  (V  notch  tip),  only  that  in  this  case  the  situation 
is  much  more  complicated  due  to  many  more  variables  involved  in  the  expression  of  the 
singular  elastic  solution  [1].  Recently,  a  general  method  for  the  accurate  and  reliable 
computation  of  all  the  information  associated  with  the  singular  solution  (generalized 
stress  intensity  foctors  -  GSlFs  -  and  associated  eigen-pairs)  in  the  vicinity  of  any  2-D 
elastic  singular  point  has  become  available  in  a  p-FEM  code  [2].  This  information  is  used 
herein  for  investigating  possible  assumptions  on  failure  criteria  from  V-notch  tips.  The 
correlation  of  a  functional  based  on  the  computed  GSIFs  and  associated  eigen-pairs  to 
e^rimental  observations  may  be  the  key  for  the  formulation  of  a  feiluie  criterion  at 
singular  points.  Herein  we  address  the  computation  of  the  required  infoimation  at 
sin^ar  points  ming  the  p-FEM,  tlw  need  for  a  functional  instead  of  a  single  value  in  the 
viciiuty  of  the  singularity,  and  preliminary  correlations  to  experimental  observations.  The 
el^c  strain  energy  density  in  tire  vicinity  of  tiie  V-notch  tip  is  suggested  as  the  foilure 
driven  functional.  Brief  explanation  is  provided.  Let  us  consider  the  vicinity  of  a  re¬ 
entrant  comer  as  shown  in  Figure  1 . 

The  elastic  strain  energy  Ufa)  density  in  Q  (the  vicinity  of  the  notch  tip)  is: 


where  [D]  is  a  dififerential  matrix,  [E]  is  the  material  matrix  and  r,  is  the  thickness  [3]. 
Based  on  the  weak  formulation,  instead  of  one  can  compute  the  linear  form  F(a) 
(involves  only  a  1-D  integration  along  the  circular  path  Fr): 
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where  f  denotes  the  traction  vector.  Let  us  denote  [f^  'u\^^dd  as  “expr-1”. 

In  the  vicinity  of  the  singular  point,  the  asymptotic  expansion  of  the  displacement  field  is: 


Were  or,, are  called  eigen-pairs  and  A,  are  the  GSIFs. 


Figure  1  -  Singular  point  vicinity:  Notations. 


Obtaining  the  stress  tensor  from  (1.3): 


- 

.  v'  'J 

(1.4) 


the  strain  dewity  is  completely  deteimined  from  the  •  JSIFs  and  associated  eigen- 
paiis  once  obtained  from  p-FEM  results:  ** 


(1.5) 

Since  in  the  vicinity  of  the  singularity,  for  R«u  the  first  term  (i  ]J  in  (1.5) 
dominates  (X,  become  larger  as  /  increase),  thus  an  excellent  approximation  for  U  is 
^ned  by  considering  the  first  term  alone.  Examining  (1.5)  it  may  be  noticed  that  the 
energy  density  is  a  measure  dependiiig  both  on  the  GSIFs  god  the  eigen-pairs, 
which  depend  not  only  on  the  loading  but  on  die  re-entrant  corner  angle  as  well. 


•c  validity  of  the  strain  energy  density  as  a  criterion  for  failure  initiation  is 
venfied  by  comparing  to  experimental  observations  reported  in  [4J. 
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p-FEM  analyses  were  performed  for  three  and  four  point  bending  re-entrant  V- 
notched  PMMA  specimens  shown  in  Figure  2,  as  well  as  sensitivity  tests  in  order  to 
establish  the  influence  of  the  load  and  support  representation  on  the  results  in  the  vicinity 
of  the  singular  point,  and  then  compared  our  results  to  the  experiments  done  by  Dunn  et 
al.  [4]. 


Figure  2  -  Four  point  and  three  point  specimens 


This  talk  will  provide  the  preliminaiy  results  of  our  analysis  and  their  conelation 
with  experimental  observations.  If  time  allows,  the  influence  of  a  rounded  V-notch  tip  on 
the  sensitiviQr  of  the  results  will  be  addressed. 
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The  recent  literature  on  variational  projection  methods  is  permeated  by  the  idea 
that  they  can  be  used  with  spatial  interpolations  not  satisfying  the  inf-sup  condition, 
often  referred  to  as  the  LBB  condition  from  Ladyzhenslmya,  Babiiska  and  Brezzi. 
In  fact,  this  IfinH  of  methods  are  often  implemented  using  equal  order  interpolations 
for  velocity  and  pressure  and  employed  as  stabilization  techniques  for  the  solution 
of  the  incompressible  Navier-Stokes  equations.  However,  severe  spatial  instabilities 
are  found  to  plague  the  numerical  solutions  obtained  by  Galerkin  finite  element 
projection  methods,  as  thoroughly  illustrated  in  [3]. 

Since  very  efficient  projection  methods  based  on  a  Galerkin— Legendre  spectral  ap¬ 
proximation  have  been  developed  in  the  last  years  [6],  it  is  worthwile  to  clarify 
whether  the  aforementioned  stability  condition  on  the  spatial  representation  of  ve¬ 
locity  and  pressure  fields  has  sensible  consequencesalso  in  the  spectral  case. 

In  the  present  paper  we  show  that  only  by  using  interpolations  of  different  order 
sudi  as  Pff-PN-2  for  velocity  and  pressure  [3]  stable  projection  method  of  a  spectral 
type  can  be  obtained  for  computing  incompressible  fiows  in  2D. 

In  particular,  the  incremental  version  of  the  projection  method  proposed  in  [2]  re¬ 
quires  to  solve  two  uncoupled  problems:  first  solve  the  viscous  step  problem 


afraid  _  -  P*'"^)  “  (u*V)u*, 


then  solve  the  incompressible  projection  step  problem  formulated  as  a  Neumann 
boundary  value  problem  for  the  pressure  increment 


(  _v2(p*+^  -  P*)  =  -(Ai)-‘  W+S 

(3)  I  a(/-^^-P*) 

[  dn  \da  ' 
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Figure  1:  Unstable  (left)  and  stable  (right)  steady  pressure  fields  for  Re  =  100 
obtained  by  an  equal  order  approximation  P^itPaio  3-  mixed  method  P40“-P38i 
respectively. 


The  elliptic  equations  for  velocity  and  pressure  have  been  recast  in  a  weak  varia¬ 
tional  form  by  the  Galerkin— Legendre  spectral  method  of  Shen  [5].  The  solution  of 
the  fully  discrete  elliptic  equations  are  calculated  by  a  direct  method  with  eigen- 
decomposition  in  both  spatial  directions.  The  nonlinear  terms  are  evaluated  by 
standard  pseudospectral  technique  based  on  Gauss-Legendre  quadrature  points. 

Figure  1  contadns  spectral  solutions  of  the  standard  driven  cavity  problem  at  Re  = 
100  obtained  by  means  of  the  non-incremental  method  with  At  =  10“^.  The  highly 
oscillatory  pressure  field  obtained  by  equal  order  spatial  discretization  Pff-Ps  is 
compared  with  the  pressure  field  computeded  by  a  compatible  mixed  approximation 
Pn-Pn-2.  The  limited  oscillations  in  the  latter,  located  near  the  boundary,  are  due 
to  the  Gibb’s  phenomenon  induced  by  the  singularity  of  the  pressure  field  in  the 
upper  comers.  These  residual  oscillations  disappear  completely  by  subtracting  the 
singular  component  of  the  solution  as  has  been  found  in  other  computations  to  be 
prtesented  in  the  final  paper. 

It  is  to  be  remarked  that  such  an  instability  has  been  observed  for  the  non-incremental 
method,  provided  the  time  step  is  sufficiently  small,  as  well  as  for  the  incremental 
method  previously  described.  Moreover,  the  same  unstable  results  have  been  ob¬ 
tained  alk)  when  the  Poisson  equation  for  the  pressure  has  been  solved  by  enforcing 
the  Neumann  boundary  condition  in  an  essential  way,  as  proposed  by  Shen  in  [5, 
page  1500]  and  in  a  private  communication,  instead  of  the  classical  weak  variational 
form. 

The  results  to  be  presented  confirm  that,  in  compliance  with  the  LBB  condition, 
the  incremental  (pressure  correction)  method  must  be  implemented  by  means  of 
interpolations  of  ffifierent  order  to  provide  a  stable  spectral  scheme  for  any  time 
step.  Such  a  method  gives  nonoscillatory  velocity  and  pressure  fields  under  ordinary 
circumstances,  without  requiring  to  introduce  any  ad  hoc  stabilization  technique. 
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The  projection  method  dates  to  more  than  thirty  years  ago  and  it  has  been  inten¬ 
sively  used  for  the  numerical  simulation  of  incompressible  viscous  flows  by  means  of 
Navier-Stokes  equations  [1,  2]. 

Notwithstanding,  a  completely  satisfactory  mathematical  theory  for  this  method  in 
its  incremental  version  has  been  exploited  only  recently  by  J.-L.  Guermond  and  L. 
Quartapelle  [3],  using  a  proper  variational  formulation. 

In  this  theory  the  projection  method  is  characterized  by  the  presence  of  two  different 
approximation  spaces  for  the  intermediate  and  final  velocity,  the  latter  belonging  to 
Such  a  velocity  can  be  eliminated  from  the  implemented  algorithm  leading 
to  a  very  simple  and  efficient  procedure,  as  shown  in  [4].  Moreover,  the  incremental 
method  is  characterized  by  a  time  splitting  error  that  allows  to  develop  a  truly 
second  order  accurate  second  order  integration  scheme. 

In  the  field  of  high  order  methods,  quite  recently,  a  new,  well  conditioned,  basis  on 
triangles  has  been  proposed  by  Dubiner  [5]  and  employed  by  Sherwin  and  Kami- 
adakis  [6,  7]  to  implement  spectral/hp  element  method  on  triangles  and  tetrahedra. 
The  discretization  has  been  applied  to  the  Navier-Stokes  equations,  employing  the 
high  order  splitting  scheme  proposed  in  [8]. 

The  aim  of  the  present  work  is  to  report  the  first  application  of  the  projection 
method  proposed  in  [4, 3],  in  the  context  of  a  spectral//ip  element  spatial  discretiza¬ 
tion  technique.  In  the  proposed  method,  differently  from  the  original  one,  the  con¬ 
vection  term  is  advanced  explicitly,  the  viscous  term  being  advanced  implicitly.  An 
efficient  method,  suitable  for  the  Direct  Numerical  Simulation  of  two-dimensional 
incompressible  viscous  flows  in  two  dimension,  is  obtained.  The  spatial  discretiza¬ 
tion  is  carried  out  by  a  mixed  element  formulation,  based  on  the  spectral//ip  element 
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method  proposed  in  [6],  employing  different  piece- wise  polynomial  spaces  for  veloc¬ 
ity  and  pressure  in  order  to  satisfy  the  inf-sup  compatibility  condition.  This  allows 
to  extend  the  Taylor-Hood  finite  elements  to  higher  degrees  as  in  [9]. 

The  stability  and  approximation  properties  of  the  proposed  method  are  assessed  by 
two  test  cases.  The  first  one,  characterized  by  a  solution  in  closed  form,  is  meant  to 
prove  the  spectral  accuracy  of  the  method  (see  Table  1).  The  second  one  exemplifies 
the  capability  of  the  proposed  method  to  cope  with  realistic  geometries. 


Degree 

Ux 

tty 

3 

2.40  10-* 

1.42  10-® 

5 

3.07  10“® 

1.90  10-» 

,  7.  , 

2,89  10-^ 

1.69  10-7 

9 

1.84  10-® 

1.01  10-® 

11 

8.32  10-“ 

4.34  10-“ 

Table  1:  L^{Q,)  error  as  a  function  of  the  polynomial  degree. 
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The  prototype  algorithm  for  solving  2I>  Poisson -and-Hdmholtz  equations  in  a  rect¬ 
angle  by  spectral  methods  is  the  diagonalization  technique  proposed  in  [5,  p.  150],  see 
also  [4,  p.  133]  and  the  reference  therein.  This  algorithm  has  been  implement^  by 
Haidvogel  and  Zang  using  Chebyshev  polynomials  and  the  tau  method  [6].  Recently, 
the  advantages  of  resorting  to  the  standard  Galerkin  formulation  in  approximating 
elliptic  problems  by  means  of  high  order  polynomials  have  been  pointed  out  [7]  and 
a  thorough  error  analysis  of  Legendre  spectral  methods  for  the  Poisson  equation  and 
the  Stokes  problem  has  been  carried  out  in  [3]. 

In  particular,  a  Legendre  basis  Po,n  bas  been  proposed  by  Jie  Shen  [8],  leading  to  a 
pentadiagonal  mass  discrete  operator  and  to  a  sti&ess  operator  equ^  to  the  iden¬ 
tity  matrix.  A  direct  Helmholtz-Dirichlet  solver  based  on  variable  separation  and 
implementing  a  diagonalization  solution  technique  has  been  derived  with  possibly 
nonhomogenous  boundary  conditions  imposed  by  an  analytical  lifting. 

The  present  authors  [1]  proposed  an  implementation  of  the  aforementioned  basis  to 
build  a  direct  solver  using  a  double  diagonalization  technique  and  a  discrete  lifting  of 
the  Dirichlet  boundary  data.  The  discrete  lifting  is  carried  out  by  suitably  augment^ 
ing  Shen’s  Po4>f  basis — a  technique  proved  successful  also  for  the  3D  problem  [2]. 
In  principle,  the  augmented  basis  should  be  able  to  accommodate  general  boundary 
conditions  including  the  pure  Neumann  problem.  Unfortunately,  in  such  an  exten¬ 
sion  the  tensor  product  structure  of  the  solver  can  not  be  preserved  owing  to  the 
particular  form  assumed  by  the  sti&ess  one-dimensional  operator  in  the  presence 
of  derivative  boundary  conditions. 

In  the  present  work,  a  Helmholtz-Neumann  direct  solver  fully  implementing  the 
tensor  product  form  is  proposed.  Following  a  suggestion  by  Shen,  we  impose  the 
Neumann  boundary  conditions  in  an  essential  way  to  develop  a  2D  spectral  solver 
based  on  variable  separation  and  double  diagonalization  technique. 

The  efficiency  and  accuracy  of  the  proposed  numerical  solver  will  be  illustrated  by 
some  numerical  tests  and  by  its  application  in  the  context  of  a  Galerkin-Legendre 
spectral  projection  method  for  the  Navier-Stokes  equations. 
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This  contribution  is  a  presentation  of  some  time  splitting  methods  applied  to  the 
incompressible  Navier-Stotes  equations.  We  describe  among  them  the  most  cur¬ 
rently  used:  Chorin-Temam  scheme  [1],  Kim  and  Moin  scheme  [2]  and  Goda  scheme 
[3]  in  their  first  and  second  order  versions.  The  special  discretization  is  carried 
out  using  the  spectral  element  method  for  the  prediction/diflPiision  and  for  the 
correction/pressure-continuity  steps.  We  discuss  the  advantages  and  drawbacks  of 
each  of  them,  with  a  particular  focuss  on  the  artificial  boundary  condition  on  the 
pressure  whidi  would  be  responsible  of  a  numerical  boundary  layer. 

We  recall  the  mathematical  results  already  proven,  we  emphasis  on  the  work  of 
Ralph  Rannacher  [4],  of  Jie  Shen  [5]  and  of  Jean  Luc  Guermond  [6].  Then,  we  show 
numerical  results  for  anal3rtical  solutions  so  that  we  can  compare  time  splitting 
algorithms  with  the  Uzawa  [7]  solution,  so  as  to  highlight  the  reliabilty  of  this  IfinH 
of  algorithms.  Moreover,  numerical  simulation  of  Navier-Stokes  flows  will  be  given 
for  different  methods  proving  again  the  efficiency  of  these  algorithms. 

Finally,  we  present  some  interesting  variants  of  each  projection  method  that  is  well 
suited  to  parallel  computations  and  provide  numerical  evidencies  of  its  efficiency  so 
as  the  substantial  gain  in  the  CPU  time. 
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This  p{q)er  is  an  overview  of  the  higher  order  functional  modeling  schemes  based  on 
finite  element,  boundary  element,  meshless,  and  wavelet  methods.  It  is  based  on  die 
research  work  of  the  principal  author  and  his  students  over  a  period  of  25  years, 
primarily,  related  to  one-D  and  two-D  problems.  The  first  discrete  numerical  mediod  is 
well  established,  so  discussions  are  mainly  focused  on  the  last  three  methods.  Apart  from 
critically  discussing  the  main  theoretical,  algorithmic,  and  implementation  characteristics 
of  the  methods,  the  p^r  expounds  on  the  relative  merits  and  demerits  of  the  methods. 
The  results  of  a  numt^  of  smooth  and  non-smooth  example  problems  are  presented. 
Finally,  the  future  potential  of  the  wavelet  method  in  discrete  numerical  mndeling  of 
continua  is  explored. 
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After  the  work  done  by  P.  Seshaiyer  and  M.  Suri,  [1]  on  the  /ip-mortar  finite  el¬ 
ements,  the  next  step  is  to  extend,  this  method  Ao  the-neariy- incompressible  elas¬ 
ticity  model  formulated  as  a  mixed  displacement-pressure  problem  so  m  to  Stokes 
equations  in  primal  velocity-pressure  variables.  Within  each  subdomain  the  local 
approximation  is  designed  on  div— stable  hp-mixed  finite  elements.  The  displace¬ 
ment  is  computed  in  a  mortared  space  while  the  pressure  is  not  subjected  to  any 
constraints  across  the  interfaces.  By  a  Boland— Nicolaides  argument  we  prove  that 
the  discrete  saddle  point  problem  satisfies  a  Babuska-Brezzi  inf-sup  condition.  The 
inf— sup  constant  is  optimal  in  the  sense  that  it  depends  only  on  the  local  (to  the 
subdomains)  characteristics  of  the  mixed  finite  elements  and,  in  particular,  it  does 
not  increase  with  the  total  number  of  the  subdomains.  The  consequences,  that  we 
are  aware  of,  of  such  an  important  result  are  twofold. 

1.  The  numerical  analysis  of  the  approximability  properties  of  the  hj^-mortar  dis¬ 
cretization  for  the  mixed  elasticity  problem  allows  to  derive  an  asymptotic  rate  of 
convergence  that  is  optimal  up  to  ^/logp  in  the  displacement. 

2.  When  the  mortar  discrete  problem  is  inverted  by  substructured  iterative  methods 
based  on  Krylov  subspaces  with  block  preconditionners,  resuming  the  proofs  of  L. 
Pavarino  and  0.  Widlund,  [2]  completed  with  those  of  Y .  Achdou,  Y.  Maday  and 
O.  Widlund,  [3]  the  condition  number  of  the  solver  depends  polylogarithmically 
on  (p,/i),  on  the  local  inf-sup  constants  and  does  not  on  the  total  number  of  the 
subdomains. 
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Application  of  displacement-based  h  elements  for  constraint  problems  of  elasticity 
and  plasticity  often  results  in  zem  or  nearly^ero-ratea  of  convergence  in  the  displace¬ 
ments,  in  energy  norm  and,  more  significantly,  in  the  stresses.  Well-known  examples 
are  the  incompressibility  constraint,  the  shear  and  membrane  constraints  appealing 
in  plate  and  shell  models  when  the  thickness  of  the  plate  or  shell  approaches  to  zero, 
and  the  constraint  resulted  by  the  boundary  layer  phenomenon  in  plate  and  shell 
models. 

Improvement  of  displacement-based  h  elements  in  linear  analysis  seems  to  have  a 
tradition;  the  applicability  of  those  numerical  ’tricks’  (reduced-selective  integration, 
hour-glass  stabilization)  in  nonlinear  analysis  is  nevertheless  questionable.  Higher- 
order  displacement-based  p  elements  are  proved  to  be  robust  in  the  displacement 
approximations,  the  computed  stresses,  however,  are  not  free  of  locking  when  prob¬ 
lems  of  incompressible  elasticity  are  considered  [1]. 

The  use  of  mixed  p  and  hp  finite  elements  is  one  of  the  possible  solutions  to  overcome 
the  aforementioned  numerical  difficulties.  A  significant  advantage  of  mixed  finite 
element  models  over  displacement  ones  is  that  they  are  robust  with  resp^  to  both 
h  and  p  extensions,  which  is  particularly  interesting  in  applications  requiring  not 
only  efficient  high-,  but  locking-free  low-order  elements  as  well  [2]. 

Although  dual-mixed  variational  principles  provide  the  possibility  of  approximating 
the  stress  space  directly  (the  variable  of  primary  interest  in  many  engineering  ap¬ 
plications),  creating  stable  finite  element  spaces  using  a  priori  symmetric  stresses  is 
not  a  simple  task  [3].  This  is  probably  one  of  the  reasons  for  an  increasing  interest 
in  dual-mixed  variational  principles  using  non-synunetric  stresses  in  finite  element 
analysis.  Development  of  higher-order  dual-mixed  p  elements  with  curved  bound¬ 
aries  seems,  however,  to  be  rather  difficult  when  the  stress  space  is  approximated 
directly  (even  if  it  is  not  assumed  to  be  a  priori  symmetric).  The  main  problem 
is  the  assurance  of  inter-element  surface  traction  continuity  for  curved  elements, 
because  the  mapping  function  between  the  master  and  actual  element  is  nonlinear. 

This  paper  presents  dual-mixed  p  and  hp  finite  element  models  for  two-dimensional 
elasticity  problems.  The  variational  formulation  is  based  on  FVaeijs  de  Veubeke’s 
dual-mixed  principle  in  terms  of  non-symmetric  stresses  and  rotations  [4].  A  priori 
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satisfaction  of  translational  equilibrium  is  achieved  by  introducing  first-order  stress 
functions  which  are  the  directly  approximated  variables  (instead  of  stresses).  This 
leads  to  a  weak  form  which  is  analogous  to  that  of  the  displacement-pressure  formu¬ 
lation  of  elasticity  and  the  velocity-pressure  formulation  of  Stokes  flow,  allowing  the 
use  of  the  stable  mixed  spaces  developed  in  [2].  Another  advantage  in  approximat¬ 
ing  first-order  stress  functions  directly  is  that  their  continuity  across  an  element 
interface  ensures  the  continuity  of  the  surface  tractions  as  well.  This  fact  makes  it 
possible  to  construct  efficient  dual-mixed  elements  with  curved  boundaries  which 
has  basic  importance  regarding  p-extensions. 

Numerical  performance  of  three  quadrilateral  dual-mixed  hp  finite  elements  is  demon¬ 
strated  through  simple  examples.  Focusing  on  nearly  incompressible  elasticity,  nu¬ 
merical  results  are  compared  to  solutions  obtained  standard  displacement-based 
hp  finite  elements.  It  will  be  shown  that  convergence  of  the  energy  norm  as  well 
as  the  stresses  is  independent  of  the  incompressibility  constraint  when  the  present 
dual-mixed  elements  are  applied,  i.e.  they  are  free  of  locking  for  both  h-  and  p- 
extensions. 
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A  number  of  physical  phenomena  are  described  by  die  Poisson  equation 

-A«  = on  Q. 

This  equation  is  a  special  case  of  the  problem  of  die  thermally  othotropic  material.  The 
two  dimensional  heat  conduction  equation  is 

^  ft  /  du  du  .  f  du  du  d\  ..  .  ^  ,*2 

(1) 

wdiere  and  ky  are  heat  conduction  coefBcients  along  the  othotropic  axes  x  and  y, 
respectively. 

A  special  case  of  the  equation  (1)  is  the  nonlinear  partial  differential  equation 


_ £_  £«" 

dx  dx  dx,  dy  dy  dy 


=  f(x,y,u)  on  Q, 


where  0  <  n  <  oo  real  number.  If  /i  =  1  the  operator  at  the  left  side  of  the  equation  (2)  is 
reduced  to  -Au. 

We  shall  consido*  the  eigenvalue  problem  of  this  type  of  equation 

^  ^u]  ^  f  ^u]  .  ,  I  j»i-l  «  .  ^ 

— - —  —  +- - —  —  +Ay  i/  =  0  m  Q,  (3) 

^xy^x  dx)  dyydy  dy) 

where  £2  is  a  convex  domain  in  R  ,  under  Dirichlet  boundary  condition 
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«1^Q=0. 

The  weak  formulation  of  the  problem  of  (3)  is 

“xV,  +  |i/y I"  ‘  UyVy^dx  =  Ain)  uvdx 
n  n 

for  any  v  e  (Q)  .  It  is  known  that  (3)  has  a  sequence  of  weak  solutions 

(A *  («),  Uk  («))  in  R  X  (Q) ,  where  0  <  A*  (n)  <  A*+i  (w) ,  ^  =  U.....  [2]. 

The  first  eigenvalue  can  be  defined  by  the  variational  principle 

j(Kr+Kr)'*' 

A,(«)=  inf  ^ - — -i - 

J|tt|  <& 

a 

The  unique  solution  (Ai{w),«|(n))  of(3)  is  called  the  first  eigei^jair  of  (3). 

We  shall  use  the  finite  element  method  to  approximate  the  first  eigenvalue  and  the 
corresponding  eigenfuncfion  to  the  nonlinear  eigenvalue  problem  of  (3)  [3].  The  finite 
element  discretization  of  the  domain  Q  is  made  by  h  and  p  version  [1].  For  the  iteration 
of  the  generalized  Rayleigh-quotient  the  Newton-j^phson  method  is  applied  such  that  the 
parameter  n  is  increased  step-by-step. 

We  have  got  ^proximate  solutions  for  the  radially  symmetric  solutions  of  the 
eigenvalue  problem  (3)  >^en  £1  is  bounded  by  the  curve 

This  curve  plays  the  same  role  in  case  of  the  nonlinear  differential  equation  (3)  as  the  unit 
circle  in  the  linear  case  vdien  n  =  l.  The  approximate  solutions  of  the  first  eigenvalues 
obtained  by  Runga-Kutta  method  will  be  compared  with  the  finite  element  solutions. 

When  n  is  bounded  by  a  square  the  analytic  solutions  are  known  for  tiie  nonlinear 
eigenvalue  problem.  The  numerical  results  obtained  by  the  finite  element  method  will  be 
compared  with  tiie  exact  values.  Illustrating  the  performance  of  tiie  h  and  p  version  we 
represent  the  relative  error  of  the  first  eigenvalue  for  different  values  of  n.  The  superiority 
of  the  /7-version  will  be  demonstrated. 

References 

[1]  Basbuska  I.,  Osbom  J.:  Eigenvalue  ProblemSy  Handbook  of  Numerical  Analysis,  Vol. 
n.,  Elsevier  Science  Publishers,  1991. 

[2] Bogn4r  G.:  Existence  theorem  for  the  eigenvalues  of  a  nonlinear  eigenvalue  problem, 
Commun.  Appl.  Nonlinear  Anal.,  4, 1997,  No.2, 93-102. 

[3]  Szab6  B.,  Babuska  I.:  Finite  Element  Analysis.  John  Wiley,  1991. 


-34- 


B-59  a 


TWO-DIMENSIONAL  INTERFACE  PROBLEMS 

Susaxme  C.  Brenner  . 

Department  of  Mathematics,  University  of  South  Carolina,  Columbia,  SC  29208, 

USA,  (brenner@math.sc.edu). 


Let  D  be  a  polygonal  domain  in  R*  which  is  divided  into  nonoverlapping  polygonal 
subdomains  . . . ,  Qj  such  that  D  = 

Consider  the  interface  problem  of  finding  u  €  Hq{Q)  such  that 

(*)  /  Pj^^-^vdx=  f  fvdx  Vi;effo(D), 

where  pi, . . . ,  pjr  are  positive  constants.  It  is  well  known  that  in  the  neighborhood  of 
a  cross  point  the  solution  u  is  not  smooth  even  if  /  is  a  (7®°  function.  In  particular, 
the  best  that  can  be  said  about  the  elliptic  regularity  of  (»>)  is  that  there  exists  a 
number  5  between  0  and  1  such  that  /  G  implies  that  u  €  n 

Depending  on  the  coeflicients  p,  and  the  angles  among  the  interfaces  at  the 
cross  points,  the  number  5  can  be  arbitrarily  close  to  0.  This  lack  of  regularity  is 
responsible  for  the  failure  of  simple  finite  element  methods  for  interface  problems. 

In  this  talk  we  will  present  optimal  order  multigrid  methods  that  take  advantage 
of  the  singular  function  representations  of  u  at  the  cross  points  and  the  extraction 
formulas  which  express  the  stress  intensity  factors  in  terms  of  u.  These  methods 
produce  approximate  solutions  of  the  form  Uh  =  T,  where  the  s^’s  are  the 

singular  functions  and  is  a  Pi  finite  element  function  on  a  (quasi)  uniform  grid 
with  mesh  size  h.  The  convergence  of  Uh  in  the  energy  norm  is  of  order  0{h)  when 
/  €  L2{CI),  and  0{h^)  when  for  I  <  j  <  J.  The  numbers  ki^  also 

approximate  the  stress  intensity  factors  at  the  rate  of  when  /  €  L2(i^)5  and 

0{h^)  when  6  for  1  <  i  <  J. 
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We  will  discuss  two  numerical  issues  in  the  electromagnetic  scattering  in  multilay¬ 
ered  media  using  mixed  potential  integral  equation  approaches.  The  first  is  the 
construction  of  high  order  RWG  basis  functions  for  curved  surfaces  and  the  second 
the  fast  calculation  algorithm  for  the  dyadic  Green’s  functions  for  multilayered  me¬ 
dia.  Basis  functions  for  combinations  of  curved  triangular  and  quadrilateral  patdies 
can  be  used.  The  Green’s  function  calculation  algorithm  is  applicable  to  arbitrary 
number  of  layers  and  frequency  range.  Several  applications  of  electromagnetic  scat¬ 
tering  in  a  multilayered  medium,  such  as  three  dimensional  discontinuities  in  VLSI 
design  and  multilayered  RF  components  and  multilayered  antennna,  are  provided. 
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Procedures  for  bounding  the  error  in  stress  in  a  finite  element  solution  have  been 
under  development  by  the  current  authors  for  over  a  decade.  Recently  the  procedures 
have  been  refined  to  allow  recovery  of  point-wise  error  bounds  on  curved  boundaries 
and  the  computations  have  become  inexpensive.  Other  researchers  have  relaxed  the 
requirranent  that  the  error  be  bounded  and  have  developed  very  successful  error 
measures  robust  enough  to  guide  adaptive  finite  element  schemes.  Reliable 
information  about  the  exact  solution  follows  after  a  sequence  of  solutions  has  been 
executed.  We  have  taken  a  more  pragmatic  approach  of  trying  to  demonstrate  that  the 
mathematics  associated  with  the  finite  element  method  is  capable  of  providing  a 
guaranteed  bound  on  the  error  in  stress  at  any  user  defined  location  in  the  solution 
domaim  Reliable  information  about  accuracy  of  the  solution  is  then  available 
following  a  single  analysis.  The  error*  referred  to  here  is  the  discretization  error 
associated  with  die  design  of  the  finite  element  mesh.  These  error  measures  can  be 
used  to  guide  refinement  but  also  to  terminate  the  analysis  after  an  initial  solution. 

Procedures  for  bounding  the  error  in  local  point-wise  derivatives  in  a  purely  Neumann 
problem  have  been  published.  Current  research  is  focusing  on  arbitrary  Dirichlet  and 
mixed  boundary  conditions  and  a  nonlinear  problem.  The  procedures  require  the 
recovery  of  a  post-processed  solution  vdiich  smoothes  the  discontinuity  in  derivatives 
across  element  boundaries.  Control  of  die  error  bounds  at  points  close  to  or  on  the 
boundary  relies  on  the  use  of  the  fundamental  solution  on  a  half  plane  for  elasticity. 
The  bounds  on  die  error  in  the  post-processed  solution  are  guaranteed. 

In  this  paper,  we  apply  the  procedures  to  linear  two-dimensional  elasticity  widi  a 
stress  concentration.  Accurate  bounds  for  the  peak  stress,  corresponding  to  a  tight 
envelope  in  which  the  exact  solution  lies,  will  be  presented.  A  new  procedure  to  solve 
patched  residual  based  problems,  to  locally  estimate  the  error,  are  being  investigated 
and  die  results  will  be  reported  at  the  conference. 
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A  method  for  connecting  a  p-element  model  to  an  h-element  model  usmg  stahc 
condensation  is  presented.  This  approach  is  based  on  equating  strain  raergies  ^  ^ 
connection  interface  of  the  two  models.  Using  least  square  fit-on  the  displacement  felds 
at  the  interface,  an  equation  is  obtained  which  can  be  used  to  transform  the  shmess 
matrix  of  the  p^lement  model  condensed  to  the  interfece  to  an  “equivalent”  shffoess 
matrix  on  the  h-side  of  the  interfece.  Similarly,  loads  acting  on  the  p-element  model  can 
be  condensed  and  transfonned  to  the  h-side* 

This  mediod  is  useful  in  cases  where  a  component  of  a  sj^tem  demands  more  detailed 
finite  element  analysis  than  the  rest  of  toe  system  to  obtain  an  accurate  solurion.  Using 
this  approach,  tois  system  can  be  broken  down  into  two  models:  toe  component  and  toe 
system  without  toe  component.  The  first  model  (the  componeiit)  is  analyzed  using  p 
elements  to  obtain  accurate  condensed  stif&iess  and  load  matrices  at  toe  connection 
intertoce.  These  condensed  matrices  can  then  be  used  as  superelements  in  subsequent  h- 
element  analyses  on  toe  second  model. 

This  method  is  implemented  as  a  new  functionality  in  an  upcoming  release  of 
Pro/Mechanica.  Numerical  examples  based  on  this  implementation  vdll  be  provided 
using  Pro/Mechanica  and  Nastran. 
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In  this  paper  we  describe  two  basic  problems  associate  to  the  introduction  of  concur¬ 
rency  in  h-refinement  for  Delaunay  meshes..  Specifically,  we  address.problems  related 
to  the  correctness  and  the  efficiency  of  parallelized  3D  adaptive  mesh  generators 
based  on  existing  scalar  guaranteed  quality  Delaunay  tetrahedralization  algorithms. 

Delaunay  triangulation  algorithms  have  been  used  very  successfully  for  unstruc¬ 
tured  grid  generation  on  sequential  machines.  Similarly,  the  Bowyer- Watson  (B-W) 
kernel  has  been  used  succcissfully  by  Chew,  Buppert  and  many  others  to  generate 
guaranteed  quality  meshes  for  2D  domains,  and  by  Shewchuk  for  piecewise  linear 
complexes  in  3D  —under  certain  constraints.  Variations  of  this  kernel  have  been 
used  by  many  others;  the  main  difference  between  the  various  algorithms  is  in  the 
treatment  of  the  domain  boundaries  (constrsunts).  IVaditionally,  the  B-W  kernel 
is  parallelized  using  domain  decomposition;  the  challenge  then  is  to  maintnit.  the 
Delaunay  property  across  the  boundaries  of  the  subdomains.  The  main  step  in  the 
B-W  kernel,  cavity  expansion,  is  based  on  a  breadth-first  search  over  the  mesh’s 
data  structures;  sequentially,  this  is  a  very  simple  task  to  accomplish.  In  parallel, 
however,  cavity  construction  becomes  much  more  difficult,  since  cavities  may  extend 
across  the  boundaries  (interfaces)  separating  adjacent  regions  of  the  mesh. 

The  expansion  of  these  multi-region  (MR)  cavities  can  be  synchronous  or  asyn¬ 
chronous,  either  halting  or  allowing  the  creation  of  new  cavities  in  regions  partici¬ 
pating  in  the  cavity  expansion.  In  an  asynchronous  implementation,  it  is  possible 
for  MR  cavities  to  intersect  (share  tetrahedra  or  faces  with)  local  cavities  or  other 
MR  cavities,  which  can  result  in  a  non-Delaunay  re-tetrahedralization  of  the  in¬ 
tersecting  cavities.  We  analyze  the  possible  intersection  cases  in  a  multithreaded 
program  execution  model  and  study  the  impact  of  concurrency  on  (1)  the  correct¬ 
ness  of  the  parallel  h-refinement  algorithm,  (2)  the  quality  of  the  triangulation,  and 
(3)  “setbacks”  in  the  progress  of  h-refinement  parallel  algorithms. 

In  addition  to  the  difficulties  of  cavities  which  extend  over  multiple  regions,  problems 
with  correctly  updating  a  region  near  by  external  boundary  in  the  face  of  concur¬ 
rency  also  arise.  For  example,  in  a  straightforward  asynchronous  parallelization 
of  Shewchuk’s  conforming  constrained  Delaunay  tetrahedralization  ^gorithm,  one 

•This  work  is  supported  by  NSF  grants:  Challenges  in  dSE,  #EIA-9726388,  Career  Award  # 
CCR-9876179,  Research  Infrastructure  #EIA-9972853,  and  IBM,  Shared  University  Researdi 
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must  do  away  with  the  ordering  on  the  splitting  of  edges,  faces,  and  elements  called 
for  by  the  original  algorithm.  Not  only  is  mesh  quality  affected  by  this  change,  but 
neither  the  correctness  of  the  resulting  mesh,  nor  the  termination  of  the  algorithm, 
can  be  proved.  These  problems  exemplify  the  need  for  the  creation  of  a  wholly  new 
parallel  algorithm  for  Delaunay  mesh  generation. 

While  the  parallelization  of  an  existing  algorithm  seems  appropriate  at  first  glance, 
our  experience  shows  that  there  is  a  need  for  a  purely  parallel  algorithm.  This  new 
algorithm  will  be  developed  keeping  in  mind  the  problems,  the  inconsistencies,  and 
the  inefficiencies  that  arise  due  to  concurrency.  The  complexity  and  the  performance 
of  our  current  parallel  implementation  clearly  demonstrate  the  need  for  such  an 
algorithm.  Our  results  also  show  that  parallel  Delaunay  meshing  based  on  the  B-W 
kernel  is  a  foreseeable  goal,  although  much  work  remains  to  be  done  to  realize  a 
usable  and  efficient  guaranteed  quality,  parallel  mesh  generator. 
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This  pq)er  presents  a  summary  of  the  authors  Ph.D.  Thesis  [1].  A  p-code  named  PHAME 
has  been  developed  at  the  LTAS  laboratory.  This  finite  element  program  is  able  to  solve 
static  and  dynamic  problems  in  linear  elasticify-bothibr  2rD  and3TD-{^>plications. 

A  posteriori  error  estimation  techniques  can  be  classified  in  three  femilies:  1)  error 
estimation  based  on  residual  of  equilibrium  and  stress  discontinuities;  2)  error  in  the 
constitutive  relation;  3)  error  estimation  based  on  gradient  recovery.  Excepting  the 
equilibrated  residual  method  [2],  all  other  methods  don’t  seem  to  be  applicable  to  high 
degree  elements.  In  most  of  the  commercial  p-codes,  there  isn’t  any  estimator 
implemented  and  the  solution  is  coiisidered  as  good  when  the  variations  of  some  user 
defined  values  between  two  h!*;rarchical  solutions  are  small  enough.  In  this  paper,  a  new 
approach  based  on  the  convergence  of  hierarchical  solutions  is  going  to  be  presented. 

In  this  work  we’ve  tried  to  take  maximal  benefit  of  the  hierarchical  formulation  to 
estimate  the  discretization  error.  The  main  idea  is  that  a  solution  is  good  when  the 
structural  response  don’t  change  too  much  when  the  discretization  is  improved.  Our  new 
approach  ba^  on  engineering  practice  is  quite  reliable.  A  measure  of  the  global  and 
lo^  quality  of  a  solution  can  be  obtained  by  comparing  two  hierarchical  solutions;  the 
error  of  p-1  solution  can  be  estimated  by  using  degree  p  solution  as  reference  solution  [4]. 
This  lea^  to  the  following  error  indicator: 

«=  /(o’,  -O’/.-lf 

a 

Where  Op  and  ap.i  are  the  degree  p  and  p-1  stress  fields,  H  the  Hooke  matrix  and  the 
volume  of  the  structure.  This  technique  gives  an  error  measure,  which  has  the  following 
characteristics: 

a)  The  error  is  always  underestimated. 

b)  The  error  repartition  is  well  estimated. 

c)  If  the  solution  is  smooth,  the  effectivity  index  (ratio  between  the  estimated  and  the 
exact  errors)  is  very  good  (close  to  1). 

d)  If  the  solution  is  singular,  the  effectivity  index  seems  to  be  independent  of  the 
mesh;  it  seems  to  depend  only  on  the  element  degree. 
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So,  if  it  is  possible  to  deteraiine  the  singularity  level,  it  is  also  possible  to  correct  the  error 
estimation.  The  problem  is  to  determine  the  singularity  level  of  the  finite  element 
solution.  Ideally,  one  likes  to  obUun  this  information  from  the  comparison  of  the  two 
hierarchical  solutions.  This  can  be  done  by  taking  benefit  from  the  following  observations 


[3 : 

a)  Singularities  occur  in  high  stress  area  and  the  singular  point  is  characterized  by  high 
stress  and  stress  gradient  values. 

is  high  in  singular  areas,  and  low  >^4iere 


b)  The  relative  stress  variation 


^ a  —G 

p  p-i 


the  solution  is  smooth. 


This  new  error  indicator  is  able  to  predict  both  the  errors  of  solution  of  degree  p  and  p-1 
and  to  deduce  the  global  and  local  convei^nce  rates  of  the  solution.  The  error  estimator 
corrected  by  an  empiric  law  has  been  validated  in  many  numerical  tests  and  has  always 
given  satisfactory  results.  Using  this“techniqne,-itiias  been  possible-to  create  an  automatic 
hp-adaptive  error  control  method,  which  is  applicable  for  2-D  and  3-D  problems. 

Many  ^plications  in  two  and  three  dimensions  will  be  presented. 
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This  is  a  progress  report  on  the  devdopmeirfr  of  the  hp  Finite  Element  Method  for 
steady-state  Maxwell’s  equations  in  both  interior  and  exterior  domains.  We  shall 
report  on  the  following  theoretical  an  implementational  issues. 

•  De  Rham  diagram  for  hp  spaces. 

•  Interpolation  error  estimates  for  the  hp  discretizations,  and  the  corresponding 
impact  on  stability  and  convergence  analysis. 

•  Convergence  analysis  for  a  two-grid  solver. 

•  Implementation  of  the  Element  Residual  Method  for  2D  hp  meshes. 

•  3Dhp90-EM  -  the  three-dimensional  hp- Adaptive  Finite  Element  Package  for 
Electromagnetics.  Examples  of  3D  simulations  (antenna  and  scattering  prob¬ 
lems). 

•  3D  implementation  of  the  infinite  element,  numerical  examples. 

•  Full  wave  analysis  of  nonhomogeneous  waveguides.  Existence  and  convergence 
results.  Examples  of  simulations. 

•  A  parallel  implementation  of  the  2D  two-grid  solver  on  Cray  3TE. 

Please  visit  my  Web  page:  (http:/ /king.ticam.utexas.edu/People/  Faculty/DEMKOWICZ/) 
for  more  information  and  copies  of  our  current  reports. 

Acknowledgement  The  work  has  been  supported  by  the  Air  Force  tmder  Grant 
F49620-98- 1-0255,  and  NSF  through  the  NPACI  project. 
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In  recent  years,  several  high-quality  software  tools  for  static  partitioning  (e.g., 
Chaco  [2],  METIS  [3],  Jostle  [6])  have  been  developed.  By  providing  convenient 
and  efficient  partitioning,  these  tools  have  enabled  the  use  of  parallel  computing 
for  a  variety  of  applications.  In  traditional  finite  element  applications,  for  example, 
static  partitioning  tools  are  used  as  preprocessors  to  divide  finite  element  meshes 
(and  computational  work  load)  among  processors. 

For  adaptive  finite  element  methods,  however,  static  partitioning  is  nOt  sufficient 
to  provide  good  parallel  performance.  The  per-processor  work  loads  in  adaptive 
methods  can  change  as  the  mesh  is  locally  refined  (/»-refinement)  or  the  degree 
of  the  approximation  is  varied  (p-refinement).  For  adaptive  methods,  dynamic  load 
balancing  is  needed  to  redistribute  work  as  a  computation  proceeds.  General  purpose 
tools  for  dynamic  load  balancing,  however,  have  not  been  readily  available.  This 
fact  is  due  to  a  number  of  differences  between  static  and  dynamic  load  balancing. 
Since  static  partitioners  are  often  run  as  preprocessors,  they  can  easily  use  file-based 
interfaces  and  data  structures  different  from  those  of  the  applications.  They  may  be 
implemented  in  serial  and  may  use  large  amounts  of  memory  and  computing  time. 
Dynamic  load-balancing  tools,  however,  must  run  along  with  the  application.  They 
must  be  implemented  in  parallel,  and  be  fast  and  memory  efficient  to  maintain  the 
scalability  of  the  application.  They  also  need  function-call  interfaces  and  methods 
to  extract  data  from  the  application’s  data  structures.  To  be  truly  general-purpose 
tools,  d3mamic  load  balancers  must  maintain  separation  between  the  data  structures 
of  the  load-balancing  algorithms  and  the  applications  using  them. 

The  Zoltan  library  [1]  addresses  many  of  the  issues  that  arise  in  implementing 
dynamic  load-balancing  algorithms.  Zoltan  is  a  general-purpose  dynamic  load¬ 
balancing  library  containing  a  suite  of  load-balancing  algorithms.  To  enable  separa¬ 
tion  of  data  structures,  Zoltan  provides  an  object-oriented  function  interface  between 
applications  and  load-balancing  algorithms.  Information  needed  by  a  load-balancing 
algorithm  is  obtained  from  the  application  through  simple  call-back  functions.  Sev¬ 
eral  different  load-balancing  algorithms  (including  interfaces  to  ParMETIS  [4]  and 
Jostle  [6])  are  available  in  the  library,  and  since  a  common  data  structure  is  not 
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required,  new  algorithms  are  easily  added  to  the  library.  Both  graph-based  and  ge¬ 
ometric  load-balancing  algorithms  are  supported.  We  are  also  incorporating  a  het¬ 
erogeneous  computing  model  into  Zoltan  to  allow  partitioning  across  heterogeneous 
sets  of  processors  or  heterogeneous  networks.  The  model  represents  a  heterogeneous 
computer  as  a  hierarchy,  with  the  entire  system’s  topology  represented  at  the  high¬ 
est  level,  subcomponents’  topology  represented  at  intermediate  levels,  and  individual 
processors’  computing  power  and  memory  capacity  represented  at  the  lowest  level. 
Our  intent  is  to  partition  both  the  parallel  computer  and  the  application  data  in 
ways  that  minimize  idle  time  on  all  processors  and  reduce  communication  across  the 
slowest  network  connections. 

In  this  talk,  we  will  describe  the  interface  and  use  of  Zoltan.  We  will  assess  its 
performance  and  overhead  costs  in  a  3D  unstructured-mesh  finite  element  code  called 
MPSalsa  [5]. 
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Many  structures  of  interest  to  aerospace  and  the  defense  sectors  have  shell  like 
geometries  that  are  stiffened  by  ribs  and/or  stringers.  The  stiffening  components 
often  vary  in  topology  and  geometry.  For  compile  three-dimensional  geometries, 
this  poses  a  great  challenge  in  the  finite  element  modeling  of  such  structures  us¬ 
ing  continuum-based  elastic  shell  elements  using  the  p- version  of  the  finite  element 
approximation.  The  issues  that  need  solution  are  two-fold: 

1.  geometric  abstraction  and  modeling  of  the  domain,  followed  by 

2.  satisfaction  of  displacement  continuity  along  the  stiffener  junctions  when  p- 
hierardiic  approximations  are  used. 

The  first  issue  is  important  because  direct  modeling  of  the  junction  details  in 
three-dimensional  space  for  individual  stiffener  configurations  is  impractical  due  to 
the  complexity  involved  in  geometric  modeling  and  generation  of  valid  curvilinear 
meshes.  Independent  specification  of  the  in-plane  and  throu^-the-thickness  inter¬ 
polation  degrees  for  the  displacement  field  lead  to  non-conforming  finite  element 
spaces  along  the  shell  junctions. 

This  presentation  will  present  modeling  scheme  where  the  shell  and  the  stiffeners 
are  geometrically  idealized  by  model  £Eu:es  and  edges,  respectively.  Each  stiffener  is 
abstracted  by  a  set  of  geometric  parameters  that  are  sufficient  to  represent  its  three- 
dimensional  geometry.  This  is  coupled  with  an  adaptation  of  the  mortar  scheme  to 
the  enforce  weak  continuity  of  the  displacement  field  along  the  stiffener  junction. 
Three-dimensional  examples  are  presented  to  show  the  use  of  the  scheme  in  solving 
problems  from  elastodynamics  and  structural  acoustics. 

tWork  funded  by  the  OfSce  of  Naval  Research. 
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The  desire  to  reduce  costs  of  composite  aircraft  structure  has  recently  increased  the 
importance  of  bonded  joint  analysis  capabilities  in  the  aeronautics  industiy.  A  new 
initiative  has  been  formed,  entitled  the  Composites  Affordability  Initiative  (CAI),  with 
the  goal  to  decrease  tiie  recurring  acquisition  costs  of  composite  aitftame  structure  by 
50%,  thereby  making  composites  more  affordable  for  tiie  next  generation  of  fighter 
aircraft.  Participants  in  the  initiative  include  the  Air  Force,  Navy,  Lockheed  Martin, 
Boeing,  and  Northrop  Grumman.  New  validated  analysis  tools  for  composite  bonded 
joints  have  been  viewed  as  the  key  enabler  for  the  ^plication  of  iimovative 
manufacturing  and  bonding  technologies  for  airframe  primary  structure.  Thus,  improved 
analysis  tools  for  design  of  bonded  composite  joints  have  become  the  focus  of  a  CAI 
analysis  tools  team.  This  paper  reviews  the  development  history  and  describes  new 
cqiabilities  of  a  p-version  Finite  Element-based  tool  developed  by  a  joint  effort  of  the 
CAI  team  and  a  software  provider. 

The  capability  of  performing  rapid  parametric  sizing  studies  was  one  of  the  central  goals 
of  tiie  team.  It  was  desired  to  develop  parameterized  solutions  for  various  bonded  joint 
configurations,  and  validate  failure  predictions  witii  CAI  test  data,  These  solutions  would 
then  be  used  by  non-Finite  Element  experts  to  perform  rapid,  preliminary  design  sizing  of 
composite  bonded  joints.  Bonded  joint  analysis  tools  previously  existed  in  the 
aeronautics  industiy  (such  as  the  industiy  standard  A4EI  code),  but  are  limited  to  simple 
joint  geometry  (typically  single,  double,  and  step-lap  joints)  due  to  the  one  dimensional 
nature  of  the  solution.  Much  more  complex  joints  are  being  utilized  to  achieve 
affordability  goals,  thus  solution  techniques  require  at  a  minimum  two-dimensional 
continuum  elasticity  analysis.  Many  finite  element  studies  utilizing  plane  strain  analysis 
appear  in  the  literature,  all  utilizing  h-element  finite  element  models.  Many  utilized 
nonlinear  models  of  the  adhesive,  and  others  included  geometric  nonlinearity.  Different 
methods  for  handling  the  so-called  **spew”  fillet  geometry,  which  is  formed  at  the 
adhesive  ftee  edge,  have  also  been  considered.  Uncertainty  of  the  shape  of  this  region 
has  also  been  studied.  Since  singularities  exist  in  the  elasticity  solution  at  free  edges  of 
composite  layer  inter&ces,  and  also  at  reentrant  comers,  consideration  of  the 
convergence  issues  and  error  control  in  these  regions  is  important  It  is  this  author’s 
opinion  that  more  attention  to  error  control  in  these  regions  is  required,  since  failure  tends 
to  initiate  in  the  vicinity  of  tiiese  singularities.  Other  important  issues  include  the  need 
for  a  modem,  user-friendly  interftice,  abilities  for  the  user  to  build-in  failure  criteria. 
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advanced  stress/strain  extraction  capabilities,  and  insensitivity  to  very  large  aspect  ratios 
caused  by  the  modeling  of  thin  composite  layers. 

The  CAI  analysis  tools  team  completed  an  extensive  study  of  the  state-of-the«art  in  the 
area  of  composite  stress  analysis  tools,  as  would  be  j^)plied  to  failure  analysis  of 
composite  bonded  joints.  The  software  StressCheck,  developed  by  Engineering  Software 
Research  and  Development,  Inc.  (ESRD),  was  identified  as  a  potential  provider,  and 
following  an  evduation  phase,  was  selected.  This  software  offers  a  p-  and  hp-element 
formulation,  which  satisfies  the  need  for  error  control,  aspect  ratio  insensitivity,  and  also 
includes  grometric  aiul  material  nonlinearity.  StressCheck  already  had  an  electronic 
handbook  interface,  which  contained  many  of  the  required  parameterization  features. 
Once  StressCheck  became  an  approved  CAI  analysis  tool  for  further  development,  CAI 
^  ESRD  worked  closely  to  define  StressCheck  modifications  required  for  the  bonded 
Joint  applications.  ESRD  has  implemented  CAI  requested  software  improvements 
throu^  both  industry  contract  fimding  and  AFOSR  sponsorship.  This  paper  will  review 
the  mt^ifications  implemented  into  StressCheck,  and  discuss  an  example  double-1^  joint 
a^ysis.  How  die  modeling  philosophy  for  this  joint  differs  firem  previous  work  will  be 
diwuss^  including  the  handling  of  singularities,  and  error  control.  Actual  Mure 
criteria  implemented,  including  stress/strain  extraction  logistics  will  not  be  reviewed,  but 
examples  of  the  general  enhanced  capabilities  will  be  discussed. 

This  p^r  also  illustrates  an  example  of  a  very  fiuitful  effort  of  the  CAI  team  in 
ideritifying  and  filling  an  analysis  tool  requirement.  An  important  conclusion  is  the  need 
for  industiy  users  to  work  closely  with  software  providers  to  precisely  tailor  software  to 
the  requirements  of  the  problem  and  the  focus  user  group.  This  experience  highlights  the 
kq^  link  being  forged  between  aerospace  industiy  researchers  and  software  providers  to 
develop  next  generation  analysis  tools. 
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We  deal  here  with  a  high-order  formulation  of  hyperbolic  conservation  laws  using  the 
discontinuous  Galerkin  method  (DGM).  Problems  that  we  intend  to  solve  are  inher¬ 
ently  transient  and  have  discontinuities  so  that  the  DGM  is  an  appealing  alternative 
to  stabilized  finite  elements  or  finite  volumes. 

We  first  propose  a  choice  of  an  orthogonal  basis  for  triangles  and  tetraedrons  in 
order  to  diagonalize  the  mass  matrix.  This  makes  explicit  time  integration  quick 
and  e^y.  This  basis  is  simply  constructed  by  applying  a  modified  Gram-Schmidt 
to  an  initial  complete  polynomial  basis. 

Monotone  methods  (i.e.  methods  that  converges  on  a  non  oscillatory  manner)  for 
solving  conservation  laws  are  of  first  order  of  accuracy.  One  can  use  piecewise  con¬ 
stant  approximations  of  the  conservative  variables  and  achieve  monotonicity.  Hence, 
this  method  will  give  poor  accuracy  in  smooth  regions  of  the  fiow  and,  moreover^ 
shocks  will  be  heavily  smeared  due  to  large  amount  of  numerical  dissipation.  We 
propose  here  a  solution  limiter  that  is  applied  after  each  time  step  so  that  oscillations 
are  removed  from  the  solution. 

For  the  adaptivity,  we  propose  an  optimization  based  methodology  for  being  able 
to  adapt  both  element  sizes  (h)  and  polynomial  degrees  (p).  We  first  defin^  what 
is  an  optimal  h  -  p  distribution  on  the  domain.  The  objective  function  is  obviously 
the  number  N  of  degrees  of  freedom  in  the  problem:  we  minimize  N  while  keeping 
the  error  under  a  givel  level. 

Finally,  we  present  numerical  results  for  several  computational  examples  in  two  and 
three  dimensions  using  parallel  tools  available  at  RPI. 
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Adaptivity  and  parallel  computation  are  essential  to  achieve  solutions  to  compu¬ 
tationally  demanding  three-dimensional  problems.  By  automatically  refining  and 
coarsening  finite  element  meshes  (h-refinement),  relocating  meshes  (r-refinement) 
and/or  varying  the  method  order  (p>refinement)  so  that  the  effort  is  concentrated 
where  resolution  is  needed,  adaptive  methods  can  increase  time  and  space  efficiency 
relative  to  traditional  techniques.  Parallel  FEM  computations  are  typically  done  by 
decomposing  the  spatial  domain  and  distributing  the  underlying  mesh  and  solution 
data  across  the  system.  With  adaptivity,  however,  meshes  and  methods  change,  and 
the  data  structures  must  support  dynamic  mesh  modification  and  migration.  Archi¬ 
tectural  variations  further  complicate  matters;  thus,  an  efficient  data  distribution  in 
a  shared-memory  environment  may  cease  to  be  so  in  a  distributed-memory,  heteroge¬ 
neous,  or  hierarchical  environment.  Partitioning  must  be  done  with  some  knowledge 
of  the  computational  environment  and  we  describe  a  hierarchical  partition  model 
that  makes  such  information  available. 

The  large  range  of  problem  sizes  to  be  addressed  and  user  computing  environments 
available  requires  software  that  runs  efficiently  on  computers  ranging  from  single 
workstations  to  the  largest  parallel  supercomputers.  As  an  example,  consider  sev¬ 
eral  common  parallel  computing  environments.  A  program  that  is  optimized  for 
a  distributed-memory  environment  with  communication  provided  by  a  fast  switch 
might  partition  the  problem  so  that  the  workload  is  approximately  equidistributed 
and  the  interfaces  between  the  subdomains  on  Afferent  processors  are  small.  A 
RtnaP  load  imbalance  might  be  tolerated  to  reduce  communication.  Tolerating  such 
a  load  imbalance  may,  however,  not  be  worthwhile  on  an  SMP  workstation,  where 
shared-memory  communication  is  fast  relative  to  a  switch.  An  Ethernet-connected 
network  of  workstations  (NOW)  is  at  the  opposite  extreme:  communication  is  very 
expensive.  Thus,  a  larger  imbalance  may  be  accepted  if  some  communication  were 
avoided.  With  the  same  program  executing  in  each  environment,  run-time  infor¬ 
mation  is  needed  to  determine  the  optimal  relationship  between  a  balanced  compu¬ 
tational  load  and  minimal  commimication.  Complexity  increases  further  with  het¬ 
erogeneous  or  hierarchical  environments.  A  system  containing  ethemet-connected 
ATM  clusters,  each  of  which  contains  a  number  of  SMP  workstations,  has  three 
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levels  of  interconnection  network.  Two  given  processes  may  communicate  by  shared 
memory,  ATM,  or  Ethernet.  There  is  no  need  for  great  expense  in  partitioning  data 
stored  on  the  same  SMP  node.  The  ATM  interface  is  relatively  slower,  but  still  fast 
enough  to  use  an  inexpensive  partitioning  and  to  maintain  a  close  balance.  Com¬ 
munication  over  Ethernet  must  be  minimized,  possibly  by  using  a  more  expensive 
partitioning  algorithm  and  by  tolerating  a  larger  load  imbalance.  Heterogeneous 
processors  (different  processor  speeds  or  memory  sizes)  are  also  important. 

The  Rensselaer  PaHition  Model  (RPM)  [4]  allows  a  hierarchical  partitioning  of  finite 
element  data  and  provides  run-time  information  about  a  parallel  computing  envi¬ 
ronment.  RPM  uses  the  mesh  structures  of  the  SCOREC  Mesh  Database  (MDB)  [1]; 
however,  many  ideas  may  apply  to  other  systems.  MDB  provides  an  object-oriented* 
hierarchical  representation  of  a  finite  element  mesh  with  three-dimensional  regio-M 
(elements)  having  links  to  their  bounding /oqes,  which  are  linked  to  their  bounding 
edges,  which  are  linked  to  their  bounding  vertices,  making  the  structures  suitable 
for  both  h  and  p-refinement  [2,  3]. 

Each  entity  in  a  distributed  mesh  is  uniquely  assigned  to  a  partition.  The  model  is 
hierarchical,  with  each  partition  assigned  to  a  specific  process  (the  process  model), 
and  each  process  assigned  to  a  machine  (the  machine  model).  Support  for  multiple 
partitions  assigned  to  a  process  allows  the  use  of  partitions  as  a  coarser  structure  for 
load  balancing  to  (i)  encourage  locality  of  reference,  (ii)  enhance  cache  performance, 
and  (m)  provide  a  means  of  performing  out-of-core  computation  by  swapping  par¬ 
titions  in  and  out  of  memory.  The  process  model  represents  the  “live”  state  of  a  job 
and  its  computational  and  communication  performance.  The  model  maps 

processes  to  the  actual  processing  nodes,  and  maps  interprocess  communication  to 
the  appropriate  network  or  interface. 

Information  provided  by  RPM  allows  a  partitioning  procedure  to  avoid  heavy  com¬ 
munication  at  interprocess  boundaries  across  slow  interfaces.  Many  partitioning  pro¬ 
cedures  minimize  the  communication  volume  of  the  total  interprocessor  boundary; 
however,  this  is  not  equivalent  to  minimizing  the  actual  cost  when  communicatioii 
speeds  vary.  Enhancements  to  existing  partitioning  algorithms,  such  as  interprocess 
boundary  smoothing  and  multilevel  partitioning,  and  the  development  of  new  al¬ 
gorithms  are  possible  with  the  availability  of  the  information  within  RPM  and  are 
being  explored. 
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In  recent  years  the  numerical  solution  of  dissipative  partial  differential  equations 
has  received  an  inflow  of  ideas  ffom  the-theo^  of  dynamieal-systems.  An  important 
contribution  is  the  use  of  approximate  inertial  manifolds,  AIMs,  in  the  so-called 
nonlinear  Galerkin  methods  [5].  One  of  the  most  often  mentioned  advantages  of 
nonlinear  Galerkin  methods  over  classical  Galerkin  discretizations  is  their  Mgher 
convergence  rates. 

In  [3],  an  inexpensive  novel  technique,  which  we  have  called  postprocessed  me¬ 
thod,  to  increase  the  accuracy  and  computational  efiiciency  of  spectral  polynomial 
approximations  was  developed.  The  postprocess  yields  the  same  accuracy  as  a 
nonlinear  Galerkin  method,  while  the  computational  cost  is  kept  nearly  the  same 
as  standard  Galerkin,  so  that  the  postprocessed  method  is  more  efficient  than  both 
standard  and  nonlinear  Galerkin  discretizations.  The  main  idea  in  the  postprocessed 
method  is  to  calculate  the  AIM  only  once,  when  the  time  integration  is  completed, 
and  add  it  to  the  Galerkin  approximation. 

In  [4]  we  introduce  a  postprocess  of  the  spectral  element  method  [2]  for  time- 
dependent  dissipative  equations.  Let  f2  be  a  two  dimensional  domain  with  smooth 
boundary.  We  consider  equations  that  can  be  written  in  the  form: 


(4) 


du{x,  t) 
dt 

u(x,  0) 


i/Au(x,  t)  -I-  Ji(u(x,  t))  -f  /(x,  t),  0<t<T,  x€Q, 

uo(x),  x€Q,  u(x,  t)  =  0,  x€  dQ. 


where  R  can  be  a  nonlinear  convective  term  such  as  R(u)  =  (u  •  V)u,  or  a  reaction 
type  term  such  as  R(ti)  =  for  ff  a  given  smooth  function. 

The  new  approach  may  be  summarize  as  follows:  Suppose  that  the  approximation 
is  wanted  at  time  T  >  0.  We  first  compute  the  Galerkin  approximation  based  on  a 
space  of  piecewise  polynomials  of  degree  N  by  numerically  integrating  the  equation 
in  the  time  interval  [0,T].  Then,  we  postprocess  by  solving  at  the  final  time  T  a 
discrete  linear  elliptic  problem  on  a  space  of  piecewise  polynomials  of  degree  M  >  N. 
The  rate  of  convergence  of  the  resulting  method  is  proved  to  be  min(iV“*“^,  M~*), 
so  that,  the  method  possesses  a  higher  rate  of  convergence  than  the  Galerkin  method 
upon  which  the  postprocess  is  applied.  Since  the  elliptic  problem  in  the  higher  order 
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poi}rnoinial  space  is  solved  only  once,  when  the  time  integration  is  completed,  the 
overcost  of  the  postprocessed  procedure  is  nearly  negligible. 

In  this  talk  we  will  show  that  the  postprocessed  method  can  be  used  as  an  a  posteriori 
error  estimator  for  evolutionary  dissipative  equations  [1].  More  precisely,  we  will 
show  that  the  error  achieved  using  the  spectral  element  method  can  be  accurately 
estimated  by  calculating  the  or  norm  of  the  difference  between  the  spectral 
element  approximation  and  the  postprocessed  approximation  that  can  be  obtained 
from  it.  The  postprocessed  method,  used  as  an  a  posteriori  error  estimator,  reveals 
itself  not  only  cheap  and  easily  computable,  but  also  able  to  give  local  and  global 
information  on  the  error  of  the  numerical  solution. 
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Recently,  a  new  type  of  time  step  integration  algorithms  using  complex  time  steps  has 
been  proposed.  The  dynamic  responses  are  evaluated  at  specific  sub-step  locations  and 
then  linearly  combined  to  yield  solutions  at  the  end  of  the  time  step.  For  linear 
problems,  tiie  algorithms  are  higher  mder  accurate,jjnconditionally.  stable  and  have 
directly  controllable  numerical  dissipation.  It  also  maintains  the  accuracy  of  the 
particular  solutions  arising  from  the  excitation.  The  complex  time  step  method  has  been 
{q>plied  to  structures  under  blast  loading  using  a  large  time  step  and  the  algorithm  can 
capture  the  high  order  loading  accurately  across  the  time  step.  The  results  are  veiy 
accurate.  However  only  linear  problems  were  considered. 

In  this  paper,  tiw  algorithms  are  extended  to  solve  nonlinear  inoblems.  As  an 
implicit  method,  the  material  state  within  the  time  interval  has  to  be  estimated  to  yield 
the  solutions  at  the  end  of  the  time  step.  To  maintain  the  dynamic  equilibrium  of  the 
problem,  nonlinear  iterative  solution  techniques  are  necessary  to  control  the  accuracy  of 
the  solutions  via  a  convergence  check.  The  nonlinear  iteration  method  chosen  is  the 
computationally  efficient  and  commonly  used  pseudo  force  method.  The  pseudo  forces 
are  Ae  internal  forces  accounting  for  the  differences  in  the  truly  nonlinear  system  and 
the  pseudo-linear  system  by  the  initial  or  fixed  state.  It  is  a  linearization  approach. 
These  forces  are  i^dated  at  each  iteration  and  treated  as  additional  excitation  on  the 
structure. 

When  the  iteration  process  is  rq)plied  to  the  complex  time  step  method,  special 
treatments  are  required  to  handle  the  pseudo  force  in  order  to  generate  accurate 
numerical  results.  The  linearization  of  the  pseudo  forces  is  removed  by  cot^airucting  the 
interpolation  function  to  reflect  the  actual  nonlinearity  of  the  system.  The :  ->^udo  forces 
are  evaluated  corresponding  to  the  computed  interpolation  fonctiotL  The  no.  linearity  is 
therefore  restored  within  the  time  step.  Since  the  complex  time  step  metix  <;i  a  sub¬ 
stepping  procedure,  the  pseudo  force  method  is  performed  independently  at  the  sub-stq> 
locations  but  the  convergence  check  is  carried  out  only  after  the  sub-step  responses  are 
combined. 

Several  numerical  examples  are  analyzed  to  demonstrate  the  applicability  of 
complex  time  step  method  on  nonlinear  problems  using  the  pseudo  force  method.  It  is 
observed  that  the  complex  time  step  method  can  be  computationally  much  more 
efficient  than  the  Newmatk  method  and  the  Runge  Kutta  method  when  very  accurate 
numerical  solutions  are  required. 
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A  fully  discrete  hp  finite  element  method  is  presented.  This  method  is  based  on  the 
Spectral  Galerkin  Element  Algorithm  which  we  describe  in  detail.  Our  method  com¬ 
bines  the  features  of  the  standard  hp  finite  element  method  (conforming  Galerkin 
Formulation,  variable  order  quadrature  schemes,  geometric  meshes,  static  conden¬ 
sation)  and  of  the  spectral  element  method  (special  shape  functions  and  spectral 
quadrature  techniques).  The  speed-up  (relative  to  standard  hp  elements)  in  the 
element  computation  results  from  our  spectral  algorithm  that  is  based  on  sum  fac¬ 
torization  and  making  use  of  properties  of  the  internal  shape  functions. 

We  present  results  on  various  problems  from  mechanics  that  show  the  computational 
advantages  of  this  algorithm  in  terms  of  CPU  time  and  convergence  rates.  We  also 
obtain  linear  speed  ups  by  performing  the  element  computations  in  parallel,  and  we 
solve  the  sparse  global  linear  S3rstem  in  parallel  that  results  from  static  condensation 
of  the  element  contributions. 
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JACOBI  SPECTRAL  METHOD 

Guo  Ben-yu 

Department  of  Mathematics,  Shanghai  Normal  University,  PR  China  200234 

(byguo@guomai.sh.cn) 


Let  A  =  {a;||a;|  <  1},  a, ^  >  -1  and  =  (l-a:)“{l+x)^.  Denote  by  (u,  v)  («.<,) 

and  llullj^ca.^)  the  inner  product  and  the  norm  of  the  space  Similarly,  we 

define  the  space  H^ia^){A)  with  the  norm  |Iv||m,x(“.«  •  For  any  even  integer  r  =  2m, 


where 


=  {v|v  is  measurable  and  ||u||r,x(«.<»),yi  <  oo} 

m-1 

=  (  E  J. 


For  any  real  r  >  0,  the  space  is  defined  by  space  interpolation.  Nect, 

for  any  non-negative  integer  /x, 

=  {v\d^v  €  =  {u|u  €  0  <  k  <  fj,} 

with  the  following  norms 

fc=0 

For  any  ^  >  0,  we  define  the  corresponding  spaces  by  space  interpolations. 

Now  let  N  be  any  positive  integer,  and  Vn  be  the  set  of  all  algebraic  polynomials  of 
degree  at  most  N.  Let  be  the  L^^a,$)—  orthogonal  projection  from 
to  Vs.  * 

Theorem  1.  Ifa  +  r>lor^-|-r>l,  then  for  any  t;  e  ^?J(o.^),,*^(A),  r  >  1  and 

where  r)  =  2/x  -  r  for  /x  >  0,  and  flr(/x,  r)  =  /x  -  r  for  /ix  <  0. 

If,  in  addition,  a  =  ^,  then  <7(^,  r)  =  2/x  —  r  —  j  for  /x  >  1,  (t(/x,  r)  =  |/x  —  r  for 
0  <  A*  ^  1»  and  <7(/x,  r)  =  /x  —  r  for  /x  <  0. 

If  the  coefificients  of  derivatives  of  different  orders  degenerate  in  different  ways,  then 
we  can  compare  the  approximate  solutions  with  the  exact  solutions  in  certain  Hilbert 
spaces.  Let  a,  ^,7,5  >  -1, 

=  {v\v  is  measurable  and  l|v||i,a,4,-y,tf  <  oo} 


-58- 


where 


For  0  <  /i  <  1,  we  define  the  space  with  the  norm  ||v||^,a,^,7,«  by  space 

interpolation. 

Let 

aa,^,7,i («» w)  =  (5x«,  +  (u,  v);j(7.<),  Vu,  V  € 

The  orthogonal  projection  Pj;r^^,y,s  is  such  a  mapping  that  for  any  v  € 

^a,P,'r,s{P}f^^,y,5‘^  ~  =  0>  €  TV- 


Theorem  2.  If  a  <  7+2  and  0  <5+2,  then  for  any  v  €  r  >  1, 

If,  in  addition,  a  <  7  + 1  and  ^  <5+1,  then  for  all  0  <  /i  <  1, 


If  the  functions  vanish  at  one  or  two  extreme  points,  we  can  define  other  orthogonal 
projections,  see  Guo  [1,  2] 

We  now  turn  to  some  applications.  We  first  consider  the  following  problem 
-dx{k{x)dxU{x))  +  b{x)U{x)  =  f{x),  xe  A, 

where  k(x)  >  0,b(x)  >  0,  and  k{x)  ~  x^**'^\x),b{x)  ~  as 

fx|  ->  1.  Then  we  can  use  the  Jacobi  spectral  method  to  solve  this  problem  numer¬ 
ically. 

Next,  let  A  =  {2/I  0  <  y  <  00}  and  consider  the  logistic  equation 

f  dtV{y,  t)  -  ^^y,  t)  =  V{y,  *)(1  -  F(y,  t)),  y  <€  A,  0  <  i  <  T, 

J  V(0,t)=J[ime-2»5,y(y,t)=0,  0  <  i  <  T, 

[  V(y,0)  =  Vo(y),  y€A. 

Let  y{x)  =  -21n(l  -  a;)  +  2  In 2,  U{x,t)  =  V{y(x),t)  and  Uo{x)  =  Vo(y(x)).  Then 

{dtU{x, t)  -  i(l  - a:)a,((l  - x)dxU{x, t))  =  U{x, i)(l  -  [/(x, t)),  xeA,0<t<T, 
U{-1,  t)  =  lim(l  -  xfdxU{x,  t)  =  0,  0  <  t  <  T, 

U{x,Q)  =  Uq{x). 

We  can  also  use  the  Jacobi  spectral  method  for  this  problem.  We  refer  to  Guo  [2] 
for  detail. 
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MULTI-P  PROCESSES  AND  PRE-CONDITIONERS 

Xian-Zhong  Guo  and  LNorman  Kats 

Dept  of  Systems  Science  and  Mathematics,  Washington  Univ.,  St  Louis,  MO  63130, 

USA  (guo@zach.wustl.edu) 

Ning  Hu 

Endocardial  Solutions  Inc.,  St  Paul,  Minnesota,  USA 


A  natural  analogy  to  the  multi-grid  method  vdiich  is  based  on  the  h-version  of  the  finite 
element  mediod  is  the  multi-p  method  which  is  based  on  the  p-version  and  hierarchic 
basis  functions.  Multi-p  processes,  including  standard  V-cycles,  m^fied  V-cycles, 
varying  V-cycIes,  and  fully  multi-p  methods,  can  be  used  as  pre-conditioners.  Methods 
for  enhancing  the  performance  of  multi-p  preconditioners  by  parallelization  or  replacing 
Gauss-Seidel  smoothers  with  damped  Jacobi  smoothers  are  described.  Results  of 
numerical  experiments  are  presented. 
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P-ADAPTIVE  TECHNOLOGY:  PERSPECTIVES  ON  ITS 
IMPACT  AND  FUTURE  GROWTH 

Marc  Halpem 

D.H.  Brown  Associates  Inc.,  222  Grace  Church  st.,  Port  Chester,  NY  10573,  USA 

(halpem@dhbrown.com) 

The  impact  of  p-adaptive  technology  reaches  &r  beyond  its  contribution  to  automatic 
detection  and  control  of  numerical  errors-in~tinite  element- models.  From  a  technical 
standpoint,  p-adaptivity  opens  opportunities  to  address  previously  intractable  challei^es 
in  diverse  disciplines  ranging  from  biomechanics  to  microelectromechanical  systems 
(MEMs).  On  a  broader  perspective,  experience  with  p-adtqitivity  has  led  to  an  increased 
role  for  CAE  simulation  in  the  mechanical  design  process.  This  improved 
understanding  of  design  requirements  will,  in  turn,  underpin  future  progress  at 
effectively  deploying  p-adaptive  capabilities. 


Over  the  past  decade,  p-adaptive  technology  has  captured  substantial  mind  share.  In  the 
late  1980s,  Noetics  emerged  as  the  first  significant  commercial  provider  of  p-ad{q}tive 
finite  element  software.  Noetics  advertised  that  its  software  would  allow  engineers  to 
focus  less  on  generation  of  acceptable  analysis  results  and  more  on  their  interpretation. 
Rasna  Corporation  then  joined  Ae  field,  identifying  and  executing  a  commercial 
strategy  to  enhance  and  promote  the  technology  through  its  Mechanica  product  line. 
Rasna  founders  believed  that  the  expertise  required  to  control  numerical  error 
represented  a  key  barrier  to  widespread  use.  By  employing  p-adaptivity  to  eliminate  that 
barrier,  they  sought  to  accelerate  tiie  growth  of  the  CAE  simulation  marketplace  fiem  a 
small  community  of  expert  technologists  to  a  much  broader  community  of  designers 
and  design  engineers.  Attracted  by  Rasna’s  rapid  growth.  Parametric  Technology 
Corporation  (PTC)  acquired  the  company  in  1 995. 


Today,  p-adaptive  Mechanica  software  represents  approximately  12%  of  mechanical 
computer-aided  engineering  marketshare  revenues  -  largely  due  to  Mechanica’s 
penetration  of  PTC’s  customer  base.  This  growth  has  compelled  other  major  suppliers 
of  finite  element  analysis  software  to  provide  some  level  of  p-adaptive  support.  All  the 
major  mechanical  engineering  analysis  software  suppliers  -  representing  more  tiian 
60%  of  the  market  -  now  support  some  level  of  p-adaptive  technology.  In  addition,  a 
1996  D.H.  Brown  Associates,  Inc.  (DHBA)  survey  suggests  that  30%  of  all  FEA 
practitioners  have  experience  with  p-adaptive  capabilities.  This  population  continues  to 
expand. 
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Despite  p-adaptive  technology’s  proven  effectiveness  in  numerical  error  control  and 
considerable  market  exposure,  its  actual  deployment  in  design  practice  has  fallen  short 
of  expectations.  More  than  just  error  control  issues,  achieving  more  extensive  use  faces 
major  barriers  such  as  the  ability  to  translate  design  questions  into  well  posed 
simulation  activities.  The  expertise  required  to  correctly  interpret  results  further 
compounds  this  challenge. 


As  a  result,  the  majority  of  designers  still  depend  heavily  on  CAE  simulation  experts 
and  design  engineers,  who  spend  at  least  half  their  time  on  engineering  analysis.  Role 
definitions  and  job  performance  metrics  for  designers  also  discourage  more  extensive 
use  of  computer-based  simulation  tools. 


Consequently,  more  effective  deployment  of  p-adaptive  tools  in  practice  will  require 
software  deployments  that  more  precisely  reflect  the  dynamics  of  design  activities  and 
capture  industry-specific  design  requirements.  The  serious  ramifications  for  design 
automation  software  based  on  these  observations  will  also  be  discussed. 
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PRECONDITIONED  ITERATIVE  SOLVERS  FOR 
COUPLED  FEM-BEM  SYSTEMS* 

Norbert  Heuer 

IWD,  Universitat  Bremen,  Germany  (heuer@iwd.uni-bremen.de) 


We  propose  and  analyze  efficient  preconditioners  for  solving  linear  systems  aris¬ 
ing  from  the  hp- version  with  quasi-uniform  meshes  of  the  finite  element/boundary 
element  coupling. 

For  our  first  strategy  we  deal  with  a  weak  formulation  that  gives  rise  to  non- 
synunetric  stifihess  matrices  whose  symmetric  parts  are  positive  definite.  These 
systems  are  solved  by  the  generalized  minimum  residual  method.  We  study  a  pre¬ 
conditioner  that  amounts  to  a  block-Jacobi  method  and  a  method  that  is  partly 
given  by  a  diagonal  scaling.  The  first  method  requires  0(log*p)  iterations  to  solve 
the  systems  up  to  a  given  accuracy,  whereas  the  second  preconditioner,  which  is 
more  easily  implemented,  requires  0(plog®p)  iterations.  Here,  p  is  the  polynomial 
degree  of  the  basis  functions. 

Our  second  strategy  is  to  consider  the  symmetric  formulation  of  the  coupling  of 
FEM  and  BEM.  Here  we  have  to  deal  with  indefinite  stiflffiess  matrices  and  we 
use  the  minimum  residual  method  as  iterative  solver.  According  to  the  structure 
of  the  Galerkin  matrix  we  study  2-  and  3-blodc  pr^onditioners  corresponding  to 
Neumann  and  Dirichlet  problems  for  the  finite  element  discretization.  In  the  case 
of  exact  inversion  of  the  blocks  we  obtain  bounded  iteration  numbers  for  the  2- 
block  Jacobi  solver  and  0(h"®/V^2)  iteration  numbers  for  the  3-block  Jacobi  solver. 
Here,  h  denotes  the  mesh  size  and  p  the  polynomial  degree.  For  the  efficient  2-block 
method  we  analyze  the  infiuence  of  various  preconditioners  which  are  based  on 
further  decomposing  the  trial  functions  into  nodal,  edge,  and  interior  functions.  By 
further  splitting  the  ansatz  space  with  respect  to  basis  functions  associated  with 
the  edges  we  obtain  a  partially  diagonal  preconditioner.  The  penultimate  method 
requires  0(log*p)  iterations  whereas  the  latter  method  needs  0(plog*p)  iterations. 

Numerical  results  are  presented  which  support  the  theory. 

*Thi8  work  is  based  on  a  collaboration  with  Ernst  P.  Stephan  and,  in  part,  with 
Maischak  (both  at  Universitat  Hannover,  Germany). 
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ON  NONLINEAR  AND  fflGHER-ORDER  PLATE 
ANALYSIS  USING  THE  P-VERSION  OF  THE  FEM 

Stefan  M.  Holzer 

Infoimationsverarbeitung  im  Konstrutiven  Ingenieurbau,  UniversitSt  Stuttgart, 
Pfaffenwaldring  7,  D-705S0  Stuttgart,  Germany  (stefaiLhoizer@po.uni-stuttgartde) 

We  discuss  an  implementation  of  the  /7-version  of  the  FEM  for  three-dimensional  elastic 
problems  under  the  restriction  that  the  geometry  of  the  structure  is  created  by  sweeping 
a  two-dimensional  domain  orthogonally  in  the  third  direction  z.  Examples  of  such 
structures  include  most  of  all  civil  engineering  structures,  like  plates,  slabs,  pillars,  and 
beams.  The  z  direction  being  perpendicular  to  the  other  two  spatial  directions, 
integrations  in  z  direction  can  be  performed  independently  and  beforehand.  Therefore, 
our  model  can  also  be  viewed  as  a  two-dimensional  /^■version  finite  element  code  for 
arbitrary  hi^er-order  hierarchical  shell  (plate)  models.  The  user  may  even  select 
different  shell  models  locally.  The  selection  of  a  specific  shell  model  corresponds  to  a 
semi-discretization  in  thickness  direction.  In  the  in-plane  discretization,  the  user  may 
also  select  an  arbitrary  combination  of  h  and  p.  We  discuss  the  influence  of  the  selection 
of  the  plate  model  on  the  solution,  requirements  for  the  in-plane  discretization, 
generalization  to  multi-layer  shell  kinematics,  curved  shells,  and  practical  iqiplications. 


In  the  second  part  of  the  talk,  we  discuss  the  application  of  the  /7-version  of  the  FEM  to 
elastoplastic  problems.  We  start  form  plane  problms  and  extend  the  discussion  to  the 
higher-order  plate  models.  Discussion  is  limited  to  J2  plasticity.  Time  permitting,  the 
talk  may  be  accompanied  by  practical  demonstrations. 
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SPLINE  VARIATIONAL  THREE  DIMENSIONAL  STRESS 
ANALYSIS  OF  LAMINATED  COMPOSITE  PLATES  WITH 

OPEN  HOLES 

E.V.  larve 

University  of  Dayton  Research  Institute,  300  College  Park,  Dayton  OH,  45469,  USA 

(Endel.Iarve@wpafb.af.mil) 


Independent  polynomial  spline  approxim^ion  of  displacement  and  interlaminar  tractions  is 
proposed  for  stress  analysis  in  laminates  with  open  holes.  Spline  approximation  eliminates 
the  inter  element  compatibility  problems  leading  to  unsatisfactory  finite  element  results  in 
the  presence  of  field  singularities.  Spline  rqjproximation  offers  continuity  of  displacement, 
strain  and  stress  fields  within  homogeneous  domains  preserving  at  the  same  time  the 
advantages  of  local  ^proximation,  such  as  sparsity  of  the  resulting  system  of  equations. 
Three  dimensional  full  field  solution  is  obtained.  Converged  stresses  and  consistent  toundaiy 
conditions,  such  as  interlaminar  traction  continuity,  are  displayed.  A  close  fprm  asymptotic 
solution,  valid  in  the  vicinity  of  the  hole  edge  at  the  interface  of  two  orthotropic  plies  of 
arbitrary  thickness  has  been  developed  to  verify  the  spline  approximation  based  full  field 
solution.  Excellent  agreement  has  been  observed  for  interlaminar  stresses  in  a  [4S/-4S]s 
AS4/3S01-6  laminate  under  uniaxial  tension.  The  polynomial  spline  ^proximation,  ideally 
suited  for  problems  concerned  witii  the  singular  solution  behavior,  has  been  applied  to  three 
dimensional  stress  analysis  in  practical  composites  containing  tens  of  plies . 


Related  Publications 

larve,  E.  V.,  “Asymptotically  Exact  Stresses  in  Laminates  ^th  Rigid  Fastener^,  J. 
Composite  Science  and  Technology,  in-press 

larve,  E.  V.  and  Pagano,  N.J.,  *^ingular  Full-field  Stresses  in  Composite  Laminates  with 
Open  Holes”,  Int.  J.  Solids  Structures,  in  press 

larve,  E.  V.,  1997,  "Three-dimensional  Stress  Analysis  in  Laminated  Composites  vtith 
Fasteners  Based  on  the  B-spline  Approximation”,  Composites  Part  A,  28A,  pp.S59-571 
larve,  E.  V.,  1996,  "Spline  Variational  Three  Dimensional  Stress  Analysis  of  Laminated 
Composite  Plates  With  Open  Holes",  Int.  J.  Solids  Structures,  44,  No.  14,  pp.  2095-21 18. 
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P-VERSION  FEA  IN  THE  AEROSPACE  INDUSTRY 
REAL  WORLD  APPLICATIONS  AND  CHALLENGES 

Tim  Jackson 

Engineering  Consultant,  Advanced  Enterprise  Solutions, 1809  East  Dyer  Road,  Suite 
313  Santa  Ana,  CA  92705,  USA  (^aclcson@aes4solutions.com) 

There  are  many  classes  of  problems  in  the  aerosp^e  industiy  for  \^ch  p-version  FEA 
is  well  suited.  This  paper  highlights  a  broad  range  of  problems  that  were  solved  using 
p-version  FEA  at  two  large  aerospace  companies.  The  method  ^ves  die  aerospace 
industiy  a  unique  and  productive  tool  to  solve  probleins  vriiere  high  element  aspect 
ratios  are  unavoidable,  easy  to  use  nonlinear  analyses  are  desired,  and  multi-material 
and  composite  analyses  involving  strong  stress  singularities  are  encountered.  Example 
problems  wll  be  shown  and  discussed  in  a  general  sense,  and  emphasis  will  be  placed 
on  the  practical  iq)plication  issues  faced  by  an  analyst  working  in  industiy  today. 


Beyond  the  technical  aspects  of  the  method,  this  paper  will  also  address  the  cultural 
issues  that  may  be  encountered  vdiile  implementing  p-version  FEA.  For  example,  the 
aerospace  industiy  is  well  entrenched  in  traditional  h^dbook  methods  and  FEA  c(^. 
Breaking  into  this  environment  with  revolutionary  analysis  tools  takes  a  great  deal  of 
patience,  persistence,  and  nurturing.  Successful  implementation  strategies  will  be 
discussed. 


While  successful  on  many  classes  of  problems,  there  are  some  problems  in  tfie 
aerospace  industiy  that  the  p-version  method  cannot  currently  solve  or  handle  with 
efficiency.  This  paper  will  discuss  anticipated  advancements  and  vriiy  they  are 
important  to  die  industiy.  Examples  of  needed  features  are:  automated  or  semi- 
automated  meshing  algorithms  for  diin  aerospace  structures,  robust  methods  to  connect 
and  analyze  assembled  parts,  and  the  ability  to  combine  and  analyze  unified  part 
domains  hairing  dissimilar  meshes. 


Beyond  the  technical  wish  list  items  is  the  need  to  continuously  educate  the  engitw^ng 
community  on  FEA.  This  effort  should  note  the  differences  in  the  underlying 
principles,  capabilities,  and  reliability  of  the  p-  and  h-version  rinite  element  methods. 
Much  too  oftm  p-version  FEA  is  misunderstood  by  analysts  and  undentfilized  by  CAD 
software  developers. 
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COMPUTATION  OF  ELECTROMAGNETIC 
SCATTERING  WITH  A  NON-CONFORMING 
DISCONTINUOUS  SPECTRAL  ELEMENT  METHOD 


David  A.  Kopriva,  Stephen  L.  Woodruff  and  M.Y.  Hussaini 
The  Florida  State  University,  Tallahassee,  FL  32306,  USA  (dak©scri.&u.edu) 


Electromagetic  scattering  problems  are  t3rpically^  cmoputed-'in  the  time>domain  by 
approximating  Maxwell’s  equations  with  low  order  finite  difference  or  finite  volume 
methods.  These  low  order  methods  have  large  phase  errors,  which  make  them 
inefficient  when  applied  to  electrically  large  objects.  High  order  finite  difference 
methods,  which  have  better  phase  properties,  are  not  practical  because  they  are 
difficult  to  apply  in  complex  geometries  or  when  material  discontinuities  are  present 
161. 

Spectral  element  methods  have  features  that  make  them  attractive  for  scattering 
problems,  when  compared  to  finite  difference  and  finite  volume  methods.  They  ran 
be  viewed  as  fiexible  extensions  of  the  spectral  collocation  method  [1],  or  as  special 
cases  oi  h  —  p  finite  element  methods.  Like  h-type  finite  element  methods,  com¬ 
plex  geometries  arising  from  scattering  bodies  or  material  interfaces  that  generate 
disconinuities  can  be  treated  easily  with  multiple  feature-fitted  elements.  Within 
each  element,  the  solution  is  approximated  by  an  orthogonal  polynomial  expansion. 
Phase  errors  can  be  controlled  either  by  decreasing  the  size  of  the  elements  or  by 
increasing  the  order  of  the  pol}rnomials  that  approximate  the  solution  [2]. 

In  this  paper  we  solve  electromagnetic  scattering  problems  by  approximating  the 
time-domain  Maxwell’s  equations  with  a  high-order  discontinuous  spectral  element 
method  (DSEM).  The  method  is  a  collocation  form  of  the  discontinuous  Galerkin 
method  for  hyperbolic  systems.  The  solution  is  approximated  by  a  tensor  product 
Legendre  expansion  and  inner  products  are  replaced  with  Gauss-Legendre  quadra¬ 
tures.  For  linear  problems  and  straight-sided  elements,  the  method  is  equivalent  to 
the  modal  form.  It  is  more  convenient  and  simpler  to  code,  however,  when  applied 
to  non-linear  problems  or  when  the  element  sides  are  curved. 

The  fact  that  the  solutions  can  be  discontinuous  at  element  faces  makes  a  DESM 
particularly  suitable  for  the  solution  of  hyperbolic  systems  in  general,  and  of  electro¬ 
magnetic  scattering  problems  in  particular.  The  method  does  not  require  continu¬ 
ous  metrics  between  elements,  which  simplifies  the  generation  of  the  mesh.  Because 
Gauss  quadratures  are  used,  special  treatment  of  element  comer  points  is  not  re¬ 
quired.  This  simplifies  programming  significantly,  increases  efficiency,  and  permits 
the  solution  of  problems  with  sharp  comers  without  additonal  code.  Finally,  the 
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discontinuous  method  permits  the  accurate  solution  of  problems  with  material  dis¬ 
continuities  at  any  approximation  order,  as  long  as  element  faces  are  aligned  with 
the  material  interfaces  [3]. 

To  increase  flexibility  of  the  DSEM,  we  introduce  a  mortar-element  implementa¬ 
tion  [5], [4].  Mortars  provide  two  benefits.  First,  they  provide  a  means  for  coupling 
element  faces  along  which  the  polynomial  orders  differ.  This  allows  the  flexibil¬ 
ity  to  choose  the  approximation  order  within  an  element  by  considering  only  local 
resolution  requirements.  Secondly,  mortars  permit  local  subdivision  of  a  mesh  by 
connecting  element  faces  that  do  not  share  a  full  side. 

The  final  paper  will  present  computations  of  some  electromagnetic  radar  cross  sec¬ 
tions  of  bodies  for  which  exact  solutions  are  known.  These  include  perfectly  con¬ 
ducting  bodies,  dielectric  bodies,  and  dielectric-coated  perfectly  conducting  bodies. 
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The  finite  element  method  is  a  widely  used  tool  for  the  solution  of  problems  in 
structural  dynamics.  Here,  the  most  popular  time  integration  algorithms  are  based 
on  Newmark’s  method,  which  represents  a  family  of  second  order  finite  difference 
formulas  that  can  be  applied  to  solve  the  semi  discrete  equations  of  motion. 

A  more  recent  development  is  the  Discontinuous  Galerkin  method,  allowing  the 
unknown  displacements  and  velocities  to  be  discontinuous  with  respect  to  time. 
The  method  is  unconditionally  stable  and  of  higher  order  accuracy.  A  two  field 
formulation  for  elastodynamics  problems  using  independent  interpolations  for  the 
displacements  and  velocities  was  introduced  by  Hughes  and  Hulbert.  The  method 
was  later  advanced  to  general  second  order  hyperbolic  systems  and  mctended  to  a 
one  field  formulation  by  the  same  authors.  A  two  field  formulation  using  a  spe¬ 
cialized  efficient  solver  for  the  resulting  linear  equation  system  was  developed  by 
Li  and  Wiberg.  Common  to  these  methods  is  the  use  of  low  order,  i.e.  constant 
or  linear,  shape  functions  for  the  time  discretization.  SchStzau  and  Schwab  ex¬ 
tended  the  Discontinuous  Galerkin  method  to  use  higher  order  shape  functions  in 
time.  They  provide  the  mathematical  framework  for  the  hp>-version  of  the  time  in¬ 
tegration  algorithm  in  the  context  of  first  order,  parabolic  systems  using  a  two  field 
formulation. 

In  our  present  contribution,  we  apply  the  idea  of  p-extensions  in  time  to  the  two 
field  formulation  for  structural  d}rnamics  and  wave  propagation  problems,  i.e.  to 
second  order,  hyperbolic  ^tems.  After  introducing  the  basic  notation,  the  high  or¬ 
der  time  Discontinuous  Galerkin  method  with  hierarchical  shape  functions  in  time  is 
outlined.  Following  ideas  of  Schotzau  and  Schwab,  we  develop  an  efficient  algorithm 
to  decouple  the  global  time  equation  system  allowing  for  a  stepwise  subsequent  so¬ 
lution.  In  numerical  examples  the  efficiency  of  the  method  is  compared  to  results 
from  Newmark’s  method  and  limitations  of  the  current  two  field  formulation  are 
discussed.  Finally  we  show  the  coupling  of  the  generically  implemented  time  inte¬ 
gration  algorithm  to  a  commercial  finite  element  package. 
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The  Finite  Element  Method  (FEM)  is  one-  of  the  primary-  tools  in  advanced  fetigue 
analysis  for  assessing  the  stresses  of  complex  machined  parts  under  cyclic  loading.  Since 
the  fatigi»  life  is  very  sensitive  to  stress  concentration,  high  accuracy  of  FEM  stress 
solutions  is  always  reqmred.  The  utmost  importance  of  any  fetigue  stress  analysis  using 
FEM  is  to  accurately  locate  the  maximum  stressed  area  and  determine  tiie  notch  stresses. 
The  part  may  be  subjected  to  complex  multi-component  loading  time  histories  and  in 
many  cases,  tiie  notch  stresses  are  multi-axial.  For  this  purpose  the  p-version  FEM  has  a 
great  advantage  over  the  h-version  FEM.  This  paper  summarizes  the  applications  of  the 
p-version  FEM  software  (MECHANICA,  STRESSCHECK)  in  tiie  notch  stress 
computations  for  crack  initiation  life  prediction  of  aircraft  structures.  The  model 
solutions  were  validated  using  a  procedure  that  includes  measurements  and  judgement. 
Ten  practical  cases  are  presented  covering  various  parts  of  the  aircraft.  Some 
comparisons  of  the  p-version  FEM  solutions  versus  solutions  from  other  methods  are  also 
presented.  Fatigue  life  computation  accuracy  due  to  variation  with  model  results  is 
presented.  Recommendations  are  presented  for  future  FEM  development 
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We  study  multiplicative  Schwarz  methods  for  the  hp-version  with  quasi-uniform 
meshes  of  the  boundary  element  method  applied  to  hypersingiilar  and  weakly  sin¬ 
gular  integral  equations  on  open  curves  and  open  surfaces. 

Proving  a  strengthenend  Cauchy  Schwarz  inequality  for  our  2-level  subspace  de¬ 
composition  and  utilizing  results  for  the  corresponding  additive  Schwarz  method 
we  show  that  the  rate  of  convergence  for  our  2-level  multiplicative  Schwarz  method 
grows  only  like  1  —  Clog”* p  towards  one  as  p  increases.  As  a  consequence,  the  num¬ 
ber  of  iterations  necessaiy  to  achieve  a  given  accuracy  of  the  discrete  solution  grows 
only  logarithmically  in  p.  Here,  p  is  the  polynomial  degree  of  the  basis  functions. 

Computational  results  are  presented  for  the  hp-version  of  the  Laplacian  and  the 
Helmholtz  operator  which  support  this  theory. 

^This  work  is  based  partly  on  a  collaboration  with  Ernst  P.  Stephan  (Universitat  Hannover, 
Germany)  and  Thanh  Ikan  (Australian  National  University,  Canberra,  Australia). 
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Introduction 

One  of  the  principal  methods  that  has  been  developed  for  error  estimation  in 
elastostadcs  relies  on  knowledge  of  both  a  conforming  or  kinematically  admissible 
solution,  and  an  equilibrating  or  statically  admissible  solution  [2,4].  TTiese  lead  to 
bounds  to  global  errors,  and  if  they  are  close  to  each  other  then  very  effective  error 
are  possible.  Whilst  dual  solutions  may  be  derived  fix>m  separate  analyses  of 
complete  models,  it  is  also  possible  to  derive  one  solution  from  the  other  using 
localised  analyses  \rith  the  aim  of  obtaining  effective  aior  estimates  at  lower 
computational  efforts.  Dual  solutions  derived  in  this  way  have  been  termed  non- 
conventional,  and  this  paper  is  principally  concerned  with  the  methods  of  deriving 
such  solutions  using  hybrid  elements  which  are  also  non-conventional  [3].  An 
alternative  approach  is  dso  proposed  based  on  a  Trefrtz  patch  recovery  method  for 
locally  smoothing  displacement  or  stress  fields,  and  for  providing  estimates  of  local 
errors  [10].  The  patches  may  consist  of  conforming  or  equilibrating  elements  such  as 
(he  hybrid  type. 

Blements  are  considered  in  the  context  of  modelling  2D,  3D,  or  plate  problems 
where  furtiier  simplifications  are  made  in  modelling  through  thickness  behaviour.  The 
primary  model  is  considered  with  a  variety  of  possible  types  of  elements  all  of  >«4iich 
have  a  /^-type  formulation  [1],  ^ere  p  refers  to  the  degree  of  displacement  or  stress 
polynomials. 

Conforming  displacement  elements  which  include  nodal  variables  at  vertices. 

Equilibrium  solutions  are  derived  by  resolution  of  nodal  forces  followed  by  recovery 
of  consistent  tractions  [6].  These  tractions  are  then  used  as  Neumann  type  boundary 
conditions  for  local  element  analyses  -  using  hybrid  equilibrium  elements  in  place  of 
^  ori^nal  elonents.  The  force  resolution  stage  does  not  provide  unique  tractions, 
various  options  are  possible  and  the  use  of  Maxwell  force  diagrams  is  presented  as  an 
engineering  solution  for  the  2D  membrane  and  plate  bending  cases.  The  3D  case  is 
i^sented  in  graph-theoretic  terms  [5]  to  represent  the  mote  complex  local 
cormectivities.  Since  the  solutions  are  not  unique,  there  is  scope  for  optimisation 
schemes  and  an  hierarchical  form  of  flexibility  method. 

Conforming  hybrid  elements  with  side  variables 
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These  elements  may  be  viewed  as  a  reformulation  of  the  conventional  conforming 
element,  but  with  side  traction  modes  and  dual  displacements  replacing  nodal 
variables.  This  implies  that  equilibrating  tractions  are  more  directly  derived,  but  the 
presence  of  spurious  static  modes  [9]  leads  to  locking  and  tractions  which  are 
indeterminate  both  fix)m  the  static  and  Mnematic  points  of  view.  This  again  implies 
that  tractions  are  not  necessarily  unique,  and  so  there  is  scope  for  local  optimisation 
which  is  discussed  in  the  paper. 

Equilibrium  hybrid  elements  with  side  variables. 

In  this  case  the  primary  model  produces  an  equilibrating  solution  for  stress  [7].  A  dual 
solution  involves  recovery  of  compatible  displacement  fields.  A  method  has  been 
proposed  [8]  based  on  element  by  element  fitting  of  displacement  fields  and  positions, 
followed  by  smoothing  of  dsplacements  at  interfaces.  An  alternative  is  proposed  in 
this  p^)er  based  on  retaining  the  interface  displacements  derived  finm  the  h3rbrid 
equilibrium  model  and  replacing  the.^original-elements.  with,  higher  degree  hybrid 
displacement  elements.  In  this  way  the  ends  of  the  interfaces  are  made  to  corifonn 
wifii  a  compatible  displacement  field,  .^ain  solutions  are  not  unique  for  the 
di^lacements  of  the  vertices,  and  there  is  scope  for  local  optimisation  of  these 
displacements.  The  question  of  spurious  kinematic  modes  is  also  ^dressed. 

Trefftz  patch  recovery  scheme  for  stresses  or  displacements. 

A  patch  recovery  scheme  based  on  Trefftz  solutions  for  smoothing  of  discontinuous 
displacement  and/or  stress  fields  has  been  proposed  [10].  This  is  discussed  in  the 
context  of  the  primary  models  considered  a^ve  with  numerical  examples  for 
Reissner-Mindlin  plates  with  Trefftz  fields  iq)  to  degree  4.  The  Treffiz  formulations 
for  a  patch  are  presented  as  for  a  “hermaphroditic”  patch  element  to  emphasize  fire 
dual  features  of  Trefftz  elements.  Local  error  estimation  in  a  patch  using  the  Trefitz 
solution  as  a  reference  solution  is  also  considered 
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For  the  mathematically  ensured,  cost-effective,  flexible  and  automatic  computa¬ 
tion  of  structural  mechanical  problems  with  error  bounds  and  tolerances,  adaptive 
finite  element  meshes  (h-adaptivity)  and  elements  with  different  Ansatz-order  (p- 
adaptivity)  as  well  as  dimension  (d-adaptivity)  are  desirable.  Isotropic  h-refinement 
and  anisotropic  p-refinement  (according  to  adequate  anisotropic  error  estimators)  is 
realized  for  the  hierarchical  Legendre  tensor  product  ansatz  functions.  Anisotropy 
is  an  important  issue  for  efficient  3D-FE-solutions  in  general.  Because  of  the  high 
numerical  effort,  the  use  of  parallel  computers  is  adequate.  This  lecture  presents 
object  oriented  data  structures  and  algorithms  which  support  these  adaptivities. 

In  most  problems  of  structural  mechanics  we  have  to  analyze  thin-walled  and  at 
least  partially  orthogonal  geometries.  These  can  be  discretized  much  better  with 
hexaeder  elements  than  with  tetraeder  elements.  We  describe  a  refinement  algorithm 
which  adapts  hexaeder  meshes  in  a  node  regular  way,  i.e.  without  hanging  nodes, 
where  degenerated  hexaeders  (-ftrieders)  appear  temporarely.  They  disappear  at 
further  refinement  which  -  together  with  adequate  numerical  integration  -  ensures 
the  full  numerical  accuracy. 

Three  topological  distinct  elements  are  needed  to  achieve  the  locking-free  transition 
from  two  dimensional  to  th  ee  dimensional  areas  of  discretizations.  Together  with 
these  hybrid  of  shell  and  3L  continuum  elements,  they  permit  the  combination  of 
both  models  in  order  to  get  the  complete  three  dimensional  solution  in  disturbed 
layers  without  the  necessity  of  using  3D-continuum  elements  for  regular  thin-walled 
subdomains.  The  transition  elements  have  continuum  nodes  with  three  degrees  of 
freedom  as  well  as  shell  nodes  with  five  degrees  of  freedom  and  couple  the  continuum 
with  the  shell  model  in  a  purely  kinematical  and  smooth  way. 

Furthermore  we  show  a  new  discretization  concept  based  on  an  hierardiically  graded 
multilevel  ansatz.  Within  this  ansatz  individual  p-orders  of  the  elements  can  be  used 
on  each  h-refinement  level.  This  is  done  in  a  decreasing  manner  from  the  coarse  to 
the  finest  mesh  level  in  order  to  balance  the  approximation  quality  of  the  total 
ansatz  and  the  computational  effort  for  solving  the  resulting  system  of  equations. 
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To  solve  the  system  of  equations  in  a  parallel  efficient  way,  iterative  multilevel  pre¬ 
conditioned  Krylov  space  methods  like  the  eg  and  Lanezos  method  in  the  symmetric 
and  the  BiCgStab  and  GmRes  method  in  the  non-synunetric  case  are  used.  With 
the  presented  hierarchical  Legendre  Ansatz  the  condition  number  compared  with 
the  Lagrangian  interpolation  polynomials  is  noticeably  smaller,  and  the  transfer  op¬ 
erator  needed,  e.g.  by  the  BPX  preconditioning,  becomes  trivial  which  improves  the 
efficiency  significantly. 

In  order  to  keep  a  general  and  simple  program  structure  for  implementing  these 
complex  algorithms  (h-,  p-,  and  d-adaptivity,  distributed  data)  they  are  coded  in 
the  object-oriented  language  C-f-H-  within  our  parallel  FE-program  ParaFep.  This 
program  also  performs  the  dynamical  load  balancing  which  is  needed  within  adaptive 
parallel  computations. 

Some  large  dimensioned  numerical  examples  illustrate  the  robustness  and  efficiency 
of  the  presented  methods. 
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Following  Irwin’s  argument,  the  energy  release  rates  are  computed  by  using  nodal 
forces  a  head  of  the  crack  front  and  relative  displacements  behind  the  crack  front. 
Furthermore,  the  quarter-point  singular  elements  were  used  for  most  finite  element 
computations.  In  this  paper,  the  energy  release  rates  are  computed  by  the  methods 
that  are  different  from  those  traditional  methods. 

Recently,  Babuska  and  Oh  introduced  a  new  approach,  called  the  method  of  auxiliary 
mapping  (MAM),  that  can  effectively  handle  singularities  in  the  framework  of  the  p- 
Version  of  the  FEM.  In  this  paper,  this  new  approach  is  modified  so  that  it  can  yield 
highly  accurate  mode-separated  energy  release  rate  for  the  delaminated  composite 
laminates.  The  results  obtained  by  this  method  are  about  10%  more  accurate  than 
most  known  results.  It  is  known  that  because  of  oscillatory  behavior  of  displacements 
near  the  crak  tip,  Gj  and  Gu  for  delaminated  composite  laminates  do  not  converge. 
However,  it  will  be  claimed  that  in  practical  point  of  view,  Gi  and  Gjt  of  interfacial 
crack  between  orthotropic  composite  materials  are  virtually  convergent.  In  this 
paper,  total  energy  release  rate  will  directly  computed  by  differentiating  the  total 
strain  energy  with  respect  to  the  crack  length. 


-78- 


0-73 


EQUILIBRATED  RESIDUAL  SPATIAL  h-  AND  p- 

error  estimators  for  elliptic,  parabolic 

AND  ELASTOPLASTIC  PROBLEMS 


Stephan  Ohnimus 

Institute  of  Structural  and  Computational  Mechanics  (IBNM),  University  of 
Hannover,  Appelstr.  9A,  30167  Hannover,  Germany 
(ohnimus@ibnm.uni-hannover.de) . 


First,  the  equilibrated  error  estimator  with  guaranteed  upper  bound  property  for 
the  error  in  the  energy  norm  will  be  given.  The  error  estimator  is  gained  by  solving 
local  Neumann  problems  with  previously  calculated  equilibrated  boundary  tractions 
regarding  to  the  current  solution.  The  equilibarted  tractions  t*  have  to  fulfill  on  the 
one  hand  the  consistency  condition,  i.e.  the  global  equilibrated  residual  7J*(v)  has 
to  be  equal  to  the  given  global  residual  72.* (v)  =  72.(v)  V  V  €  V  C  H^.  This  leads 
to  the  C^-  continuity  condition  of  the  equilibrated  boundary  tractions  in  normal 
'  direction  at  the  finite  element  boundaries.  On  the  other  hand  the  solvability  of 
the  local  Neumann  problem  leads  to  the  equilibrium  condition  for  the  equilibrated 
residuals  72^  (v)  =  0  V  v  €  Z(fte).  Herein  is  Z(f2e)  the  kernel  (rigid  body 
displacement  modes)  of  the  local  element  bilinear  form.  Examples  will  explain  the 
efficiency  of  the  error  estimator,  see  e.g.  [1,  2,  3]. 

Second,  for  parabolic  problems  a  special  norm  in  space  and  time  will  be  presented  in 
the  frame  of  semi-discrete  finite  element  method,  discrete  in  space  and  sufficiently 
exact  (quasi  continuous)  in  time.  Referring  to  this  norm  a  guaranteed  upper  error 
estimator  will  be  presented  which  used  the  results  form  the  previous  explained  error 
estimator  for  elliptic  problems,  see  [4].  This  error  estimator  is  computed  again 
by  solving  local  Neumann  problems  with  integration  over  the  time  to  take  into 
consideration  the  acciunulation  of  the  error  in  time  direction.  Two  examples  will  be 
given  for  a  ID-  and  a  2D-  parabolic  problems  in  application  to  elasticity. 

Third,  for  problems  with  elasto-plastic  material  behavior  a  guaranteed  spatial  error 
estimator  will  be  presented  under  the  assumption  that  the  ciurent  tangent  is  suffi¬ 
ciently  exact,  see  [5].  An  exapmle  will  present  the  effciency  of  the  error  estimator. 
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This  lecture  deals  with  mechanical  contact  problems  using  displacement  based 
/^■versional  finite  element  method.  In  die  contact  region  both  the  normal  stress  and  the 
stresses  may  have  singularities  as  well  as  jumps  in  their  derivatives.  The 
boundary  of  the  contact  zone  is  unknown  qpriori.  Therefore  the  locations  of  the 
singularities  in  the  normal  and  tangential  stresses  and  jumps  of  their  derivatives  are  also 
unknowiL 

Concerning  the  finite  element  discretization  we  have  problem  of  category  C  [1] 
that  is  the  in  2D  caimot  be  constructed  so  that  the  points  where  the  solution  is  not 
analytic  are  not  in  vertices.  When  the  p-version  is  used  then  the  accuracy  is  typically  high 
enough  for  the  singularities  to  induce  oscillations  in  the  numerical  solutions.  For  the 
treatment  of  axially  symmetric  contact  problems  an  adaptive  finite  element  method  has 
been  developed  by  which  the  category  C  problem  is  converted  into  category  B. 

The  mesh  is  adjusted  in  the  iteration  process  such  that  the  boundary  of  the  contact 
zone  and  the  boundary  of  the  stick-slip  zones  are  nodal  points,  allowing  the  jumps  in  the 
derivatives  to  be  represented  in  the  discretized  problem  [2].  The  accuracy  can  be 
increased  by  the  eruichment  of  the  finite  element  ^)ace  with  special  singular  fimctions 

[3]. 

In  the  present  work  it  will  be  assumed  that  tire  displacements  and  deformations 
ate  small,  the  material  of  the  contacting  bodies  are  elastic.  The  augmented  Lagrangian 
technique  is  used  for  the  solution  of  the  contact  problem. 

The  p-version  of  the  finite  element  method  allows  relatively  coarse  rnesh.  The 
positioning  of  the  nodal  points  in  the  border  of  the  contact  zone  is  performed  in  one  or 
two  phases  HapanHing  on  the  predefined  accuracy  to  be  achieved.  The  first  phase  is  rough 
positioning  of  tile  contact  points  its  aim  is  to  achieve  the  situation  when  each  of  the 
integration  points  of  the  contacting  elements  is  in  contact  In  the  second  phase  a  fine 
adjustment  of  the  location  of  nodal  points  is  performed  on  the  bases  of  monitoring 
different  indicators  (the  potential  energy,  the  smootimess  of  the  contact  jnessure/  normal 
displacement  and  contact  pressure  at  the  border  point  of  the  contact  region). 


This  woric  was  partly  supported  by  AKP-97-146,2,3/31  and  OTKA  T025172. 
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In  case  of  fnctional  contact  problem  a  nodal  point  in  the  contact  region  is 
positioned  to  the  border  point  of  the  stick  and  slip  zones  in  order  to  satisfy  the  Coulomb 
friction  law. 

The  efficiency  of  the  proposed  method  will  be  demonstrated  through  couple  of 
examples.  There  will  be  examples  in  ^Mdiich  the  light/lefr  side  or  both  sides  of  the  contact 
zone  will  be  positioned. 

Two  types  of  mechanical  contact  optimization  problems  vdll  be  investigated: 

1 .  optimization  by  the  controlling  of  the  contact  pressure  distribution, 

2.  maximization  of  the  global  mechanical  parameter  concerning  a  machine  tool 
operation. 

The  contact  pressure  generally  may  have  singularify  in  the  contact  zone  without 
shape  optimization.  Our  aim  is  to  control  the  contact  pressure  to  have  smooth  contact 
pressure  as  defined  by  a  control  function. 

The  control  can  be  described  by  the  following  formula 


X(S)  =  V(S)P^-P(S)=^0,  S€Cl^ 

>^ere  s  is  the  curve  coordinate  along  control  zone,  v(s)  -  controlling  function 
prescribed  by  the  user,  -  maximum  of  the  contact  pressure,  p(s)  is  the  contact 
pressure.  On  the  complement  subregion  of  that  is  on  die  no  control  part 
have  an  inequality  ;jf(x)>0.  We  have  prescribed  die  controlling 

function  to  be  C'  class  function,  actually  constant  in  the  middle  of  the  contact  zone  and 
Hermite  polynomials  at  the  border  zones. 

The  optimization  problem  is  a  restricted  linear  programming  problem.  For  the 
solution  a  special  iterational  algorithm  has  been  developed  vs^ch  proved  to  be  faster  than 
standard  mathematical  program  package. 

The  following  problems  will  be  analyzed: 

1 .  Contact  optimization  problems  of  cylindrical  bodies  in  overlying  case. 
Computation  of  the  stress  state  due  to  assembly  process,  considering  the  heat 
effects. 

2.  Contact  shsqM  optimization  of  cylindrical  clamping  bodies  when 
maximization  of  the  torsion  moment  transmission  is  to  be  achieved. 
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Successful  large  scale  (capable  of  solving  problems  with  millions  of  unknowns)  par¬ 
allel  adaptive  hp  codes  require  efficient  data  management  schemes,  good  adaptive 
strategies  and  fast/robust  solvers  that  can  reliably  solve  inregularly  sparse  and  often 
poorly  conditioned  linear  systems.  We  begin  by  describing  recent  work  on  develop¬ 
ing  a  code  infrastructure  that  performs  the  necessary  data  management,  partition¬ 
ing,  parallel  refinement/enrichment  and  dynamic  load  balancing.  Central  to  this 
development  is  the  use  of  a  self-organizing  data  management  scheme  based  on  bal¬ 
anced  tree  data  data  structures  on  each  processor  composed  into  a  distributed  list 
structure.  Integrated  with  this  structure  we  have  also  developed  a  class  of  solvers 
described  next. 

In  past  work,  on  these  methods  we  have  developed  iterative  substructuring  based 
solvers  and  coarse  grid  hierarchical  preconditioners  that  are  quite  efficient  for  2D 
elliptic  problems.  We  note  that  our  use  of  automatic  partitioning  algorithms  and 
unstructured  grids  precludes  the  nested  coarse  grids  that  have  been  traditionaly 
used.  Our  coarse  grid  functions  are  derived  naturally  as  linear  combinations  of  fine 
grid  functions.  For  such  problems,  we  have  proved  that  they  have  condition  numbers 
that  grow  at  rates  that  are  polylogarithmic  in  the  mesh  discretization  parameters. 
These  solvers  were  also  naturally  extended  to  the  Stokes  equations.  However,  on 
more  complex  3D  problems,  such  efficiencies  are  not  always  obtained.  In  related 
work  we  have  also  developed  a  multi-level  substructuring  type  solvers  based  on  a 
good  fill-reducing  ordering  using  a  space-filling  curve  passing  trough  the  element 
centroids.  This  ordering  also  constitutes  the  basis  for  our  data  indexing  and  data 
management  schemes.  This  solver  is  very  robust,  but  not  as  efficient  as  the  iterative 
substructuring  solvers. 

We  report  here  the  development  of  hybrid  solvers  that  combine  the  advantages  of 
both  by  automatically  switching  from  one  algorithm  to  the  other,  depending  on 
the  problem,  to  produce  robust  solvers  that  work  reliably  on  a  large  class  of  3D 
?:>roblems.  The  solver  starts  with  a  substructuring  on  the  lowest  levels  (bubble 
and  edge  function  degrees  of  freedom)  that  are  local  to  a  processor,  and  switches 
to  a  preconditioned  iterative  solver  at  a  level  of  the  hierarchy  that  require  access 
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to  non-local  data.  The  preconditioner  is  a  very  coarse  grid  and  diagonal  blocks 
corresponding  to  the  remaining  degrees  of  freedom.  If  the  solver  exhibits  difficulty  in 
converging,  we  expand  the  definition  of  the  preconditioner  to  include  more  and  more 
degrees  of  freedom  until  we  reduce  to  the  multi-level  substructuring  solver  described 
earlier.  Further,  since  the  choice  of  which  solver  is  more  efficient  on  a  multi-processor 
architecture  is  dependent  on  the  relative  costs  of  computation  and  communicaition 
it  is  possible  to  tune  these  solvers  to  particular  machine  architectures  by  selecting 
the  number  of  levels  of  substructuring  to  use. 
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This  paper  details  APES’s  experience  ^plying  the  higher  order  (^version)  finite  element 
software  StressCheck  (ESRD,  Inc.,  St  Louis,  MO,  USA)  to  highly  specialized  and 
difficult  problems  in  structural  life  predictions,  namely  in  predictions  of  iurcraft  structural 
fiitigue  life  using  crack  growth  analyses.  A  typical  crack  growth  simulation  uses  the 
geometry  and  loads  in  a  structure  to  compute  Mode  I  Stress  Intensity  Factors  (SIFs)  Ki. 
The  SIFs  are  used  to  derive  crack  growth  rates,  da/dN,  which  can  be  interpreted  as  an 
intrinsic  material  behavior,  in  addition  to  properties  such  as  Young’s  modulus  and 
Poisson’s  ratio,  for  a  broad  range  of  environments.  The  crack  growth  rate  is  used  to 
calculate  an  increment  in  the  crack  length  (da)  for  a  given  number  of  fatigue  cycles  (dN). 
The  increased  crack  lengdi  is  computed,  the  new  SIF  estimated,  and  the  simulation 
cnntiniifts  until  fiacture  of  die  structure  is  predicted.  Typically,  a  crack  growdi  simulation 
will  require  thousands  of  cycles  (dN).  Fortunately,  accurate  and  reliable  stresses  and 
stress  intensity  ftictors  derived  fiom  StressCheck  simulations,  obtained  quickly  and 
efficiently,  enable  the  basic  procedure  to  be  successfully  applied  to  development  of 
advanced  life  prediction  methods. 

Many  FEA  software  packages  could  have  been  used  to  obtain  engineering  data 
such  as  SIFs  and  stresses.  However,  the  advanced  numerical  methods  incorporated  into 
higher-order  finite  element  methods  such  as  StressCheck  allow  the  analyst  to  obtain  more 
accurate  and  reliable  results  with  convenient  error  checking,  numerical  convergence  of 
enginftpring  data  such  as  maximum  stress  and  strain  and  SIFs,  stress  contours,  and 
deformations.  In  addition,  greatly  shortened  analysis  schedules  are  realized  through  the 
use  of  procedures  such  as  the  automatic  extraction  of  SIF  parameters  by  using  the 
Contour  Integral  Method  (CIM).  The  versatility  of  higher-order  methods  is  demonstrated 
by  some  examples  that  we  have  encoimtered  in  engineering  practice — multi-site  damage 
scenarios  and  severely  corroded  plates. 

Corrosion  is  a  difficult  structural  problem  to  solve,  requiring  innovative  solutions 
in  order  to  overcome  substantial  challenges.  Often  the  engineer  will  have  a  myriad  of 
structural  effects  to  understand,  such  as  stress  risers  caused  by  corrosion  topogrsqrhies, 
bending  stresses  induced  by  corrosion  by-product  build-up  between  adjacent  parts,  and 
multiple  cracks  in  the  structure  induced  hy  corrosion  and  ftdgue  cycling.  The  advances  in 
i^•version  FEM  have  provided  the  necessary  inputs  to  crack  growth  analyses  used  to 
predict  the  effects  of  corrosion  on  aircraft  structural  life.  Applications  of  advancements 
made  in  StressCheck  have  enabled  incorporation  of  the  key  analysis  characteristics 
necessary  to  model  tiie  effects  of  corrosion  on  structural  life  predictions. 
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This  paper  presents  the  design  and  evaluation  of  an  adaptive,  system  sensitive 
distribution/load-balancing  framework  for  distributed  adaptive  grid  hierarchies  that 
underlie  parallel  adaptive  mesh-refinement  (AMR)  techniques  for  the  solution  of 
partial-differential  equations.  The  framework  uses  application  and  system  state  in¬ 
formation  to  select  the  appropriate  distribution  scheme  at  run-time.  The  flection 
is  driven  by  an  application-centric  performance  characterization  of  dyn^ic  parti¬ 
tioning  and  load-balancing  tediniques,  and  is  governed  by  rules  defined  in  a  policy 
database.  The  primary  motivation  for  the  framework  is  the  design  and  development 
of  policy  driven  tools  for  automated  configuration  and  run-time  management  of 
tributed  adaptive  applications  on  dynamic  and  heterogeneous  networked  computing 
environments. 

Dynamically  adaptive  methods  for  the  solution  of  partial  differential  equations  that 
employ  locally  optimal  approximations  can  yield  highly  advantageous  ratios  for 
cost/accuracy  when  compared  to  methods  based  upon  static  uniform  approxima¬ 
tions.  These  techniques  seek  to  improve  the  accuracy  of  the  solution  by  dynanaically 
refining  the  computational  grid  in  regions  of  high  local  solution  error.  Distributed 
implementations  of  these  adaptive  methods  offer  the  potential  for  the  accurate  solu¬ 
tion  of  realistic  models  of  important  physical  systems.  These  implementation  h<w- 
ever,  lead  to  interesting  challenges  in  dynamic  resource  allocation,  data-distribution 
and  load  balancing,  communications  and  coordination,  and  resource  management. 
The  overall  efficiency  of  the  algorithms  is  limited  by  the  ability  to  partition  the 
underlying  data-structures  at  run-time  so  as  to  expose  all  inherent  parallelism,  min¬ 
imize  communication/synchronization  overheads,  and  balance  load.  A  critical  re¬ 
quirement  while  partitioning  adaptive  grid  hierarchies  is  the  maintenance  of  logical 
locality,  both  across  different  levels  of  the  hierarchy  under  expansion  and  contrac¬ 
tion  of  the  adaptive  grid  structure,  and  within  partitions  of  grids  at  aU  levels  when 
they  are  decomposed  and  mapped  across  processors.  The  former  enables  efficient 
computational  access  to  the  grids  while  the  latter  minimizes  the  total  conunumcar 
tion  and  synchronization  overheads.  Furthermore  application  adaptivity  results  in 
application  grids  being  created,  moved  and  deleted  on-the-fiy,  making  it  is  necessary 
to  efficiently  re-partition  the  hierarchy  so  that  it  continues  to  meet  these  goals. 
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Moving  these  applications  to  dynamic  and  heterogeneous  networked  computing  envi¬ 
ronments  introduces  a  new  level  of  complexity.  These  environments  require  selecting 
and  configuring  application  components  based  on  available  resources.  However,  the 
complexity  and  heterogeneity  of  the  environment  make  selection  of  a  ’’best”  match 
between  system  resources,  application  algorithms,  problem  decompositions,  map¬ 
pings  and  load  distributions,  conununication  mechanisms,  etc.,  non-trivial.  System 
dynamics  coupled  with  application  adaptivity  makes  application  configuration  and 
run-time  management  a  significant  challenge.  Clearly  there  is  a  need  for  "smart” 
tools  that  can  automate  the  configuration  and  management  process. 

This  paper  first  presents  an  application-centric  characterization  of  distribution  mech¬ 
anism  for  AMR  applications  on  heterogeneous  (and  dynamic)  cluster  computing 
environment.  It  then  describes  the  design,  implementation  and  evaluation  of  an 
automated  application  configuration  and  management  framework  that  dsmamically 
adapts  distribution  and  load-balancing  mechanisms  based  on  current  application 
and  system  state. 
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It  is  today  widely  accepted,  that  the  p-  and  hp-version  of  the  fimte  element  method  is 
superior  to  the  classical  h-version  with  respect  to  accuraQf  and-computational  efficiency, 
if  linear,  elliptic  problems  in  two  dimensions  are  considered.  Yet,  many  researchers  still 
doubt,  that  these  advantages  can  also  be  observed  for  more  complex  non-linear  or  three- 
dimensional  problems.  It  will  be  show  in  this  contribution,  that  the  superiority  of  the  p- 
version  may  be  even  more  dramatic  in  these  cases. 

We  will  first  focus  on  the  aspect  of  geometric  modelling  in  three  dimensions.  The 
well-known  blending  function  technique  for  msqiping  the  geometry  of  p-elements  offers 
the  possibility  to  completely  separate  all  geometric  computation  involved  in  a  finite 
element  analysis  fiom  the  non-geometric  parts.  This  separation  allows  to  design  a 
distributed  software  system^  where  the  geometric  model  of  a  CAD-program  is  directly 
linked  in  a  client-server  software  structw-e  to  a  finite  element  kernel.  This  concept  offers 
the  advantage  of  using  all  state-of-the-art  CAD-techmques  like  geometric  editing  or 
parametric  design  in  a  finite  element  analysis.  The  increase  of  efficiency  for  practical 
woric  can  be  dramatic,  as  such  a  system  for  computer  integrated  engineeiing  relieves  a 
user  fiom  the  burden  to  transfer  geometric  data  fiom  CAD  to  FEA,  which  is  usually  very 
time-r.f>n«giimingj  even  if  Only  some  geometric  parameters  of  the  model  change. 

The  ptqier  will  describe  the  software  structure  of  our  integrated  system  in  detail. 
Some  basic  concepts  of  element  matrix  computation  will  be  reviewed  in  the  light  of  the 
blending  function  technique.  We  will  outline  the  advantages  of  this  approach  for  design 
processes  in  civil  engineering  ^iplications.  Special  enqihasis  will  be  laid  on  spatial 
constructions  like  shells  and  combmations  of  shell  and  three-dimensional  solid  models. 

The  second  part  of  the  psper  will  address  physically  non-linear  problems.  As  a 
model  problem  we  will  consider  the  rate  independent  plasticity  with  non-linear  isotropic 
hardening,  also  known  as  the  flow  theoiy  of  Prandtl-Reuss.  Small  strains  will  be  assumed 
and  a  vow  Mises  yield  criterion  will  be  ^plied.  As  an  integration  algorithm  for  the  rate- 
independent  equations  of  plasticity  we  will  consider  an  implicit  Euler  schane  by  applying 
a  retum-m^ping  algorithm. 

Finally,  we  will  compare  our  p-version  approach  to  the  classical  h-version  for 
some  benchmark  problems  as  well  as  for  more  complex  real  life  problems.  Applications 
for  civil  engineering  constructions  as  well  as  for  bio-mechanical  problems  will  be 
presented. 
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CURRENT  USAGE  OF  TWO  P-VERSION  FINITE 
ELEMENT  ANALYSIS  CODES  AND  PROPOSED 
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Paul  A.  Reid 
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Louis,  MO,  USA  (paul.a.reid@boeing.com) 


The  purpose  of  this  paper  is  to  present  die  current  us^e  of  Finite  Element  Analysis 
(FEA)  tools  at  Boeing  -  St  Louis  as  well  as  recommended  enhancements  for  increased 
cqiabiiity.  Two  types  of  p>version  FEA  codes  ivill  be  discussed;  StressChwk  by 
Engineering  Software  Research  &  Development  and  Pro/MECANIC A  by  Parametric 
Technology  Corporation. 

The  pi^>er  will  focus  on  conventional  strength  and  fracture  mechanics 
^plications  that  are  performed  by  structural  engineers .  Current  applications  include 
design  trade  studies,  manu&cturipg  support,  and  advanced  design  efforts.  Included  is  a 
comparative  discussion  of  strengths  and  weaknesses  of  the  various  features  available  in 
each  code. 

The  paper  will  also  propose  concepts  to  promote  future  growth  and  increased 
capability.  Discussion  is  provided  on  recommended  feature  modifications  and  topics 
for  new  analytical  tool  ctyiability.  The  purpose  of  this  proposal  is  to  initiate  dialogue 
witiiin  the  FEA  software  community  to  encourage  development  growth.  A  path  will  be 
suggested  for  both  software  designer  and  customer  to  remain  ahead  of  their  respective 
competition  by  developing  new  ways  to  increase  capability,  improve  accuracy,  and 
reduce  analysis  time. 
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A  POSTERIORI  ERROR  ESTIMATES  FOR  A 
DISCONTINUOUS  GALERKIN  METHOD  APPLIED  TO 

ELLIPTIC  PROBLEMS 

Beatrice  Rivifere  and  Mary  F.  Wheeler 
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78712,  USA  (riviere@brahma.ticam.utexas.edu) 


A  posteriori  error  estimates  for  locally  mass  conservative  methods  for  subsurface  flow 
are  presented.  These  methods  are  based  on  discontinuous  approximation  spaces  and 
referred  as  Discontinuous  Galerkin  methods.  One  can  show  an  a  priori  exponential 
rate  of  convergence  for  elliptic  problems.  Here,  the  hp  adaptivity  is  investigated  for 
flow  problems  in  2D.  Two  types  of  global  error  estimators  are  derived:  an  explicit 
L*  error  estimator  that  can  be  used  as  an  error  indicator  and  an  implicit  error 
estimator  that  is  based  on  the  norm  of  residual. 
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GERMANY 


The  first  part  of  the  lecture  is  concerned  with  a  short  review  on  regularity  results 
for  solutions  of  mixed  boundary  value  problems  for  systems  of  stationary  quasilinear 
and  semilinear  partial  dilSerential  equations  in  domains  with  conical  points,  edges 
and  vertices  [1].  We  introduce  weighted  Sobolev  spaces  with  attached  asymptotics 
generated  by  the  a^mptotical  expansions  of  solutions  of  corresponding  linearized 
problems  near  the  boundary  singularities.  We  formulate  conditions  which  guarantee 
that  the  solutions  of  the  nonlinear  problem  have  the  same  as}rmptotical  structure  as 
the  solutions  of  the  linearized  problem. 

In  a  second  part  we  apply  these  abstract  investigations  to  general  boundary  value 
problems  for  the  Navier-Stokes  equation  which  model  complicated  flow  conditions 
on  the  boundary.  The  resulting  regularity  results  are  crucial  for  an  error  analysis  of 
finite  element  methods.  Error  estimates  for  uniform  and  graded  meshes  are  derived 
[2]. 

References 

[1]  M.  Bodmiak,  A.M.  Sandig,  Local  solvability  and  regularity  results  for  a  doss 
of  semilinear  elliptic  problems  in  nonsmooth  domains,  Mathematica  Bohemica 
2-3,  124  (1999)  pp.  245-254 

[2]  M.  Orlt,  Regularity  and  error  estimates  for  general  boundary  value  problems  to 
Navier-Stokes  equations,  Thesis,  University  Rostock  (1997) 
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Jonathan  B.  Ransom 

NASA  Langley  Research  Center  Hampton,  Virginia 
Mohammad  A.  Aminpour 

Applied  Research  Associates,  Inc.Hampton,  Virginia  W.  Jeiferson  Stroud 
NASA  Langley  Research  Center  Hampton,  Virginia 

When  performing  global/local  analysis,  the  issue  of  connecting  dissimilar  meshes  often 
arises,  especially  vsdien  refinement  is  performed.  One  method  of  connecting  these 
dissimilar  meshes  is  to  use  interfiice  elements.  Curve  interftice  elements,  to  connect 
dissimilar  p-element  edges  along  a  curve;  surface  interfece  elements,  to  connect 
Hisftiitiilar  p-element  faces  over  a  sur&ce;  and  shell-to-solid  transition  interface  elements, 
to  connect  dissimilar  p-element  edges  with  p-element  ftices,  have  been  implemented  in 
MSC/NASTRAN.  The  background,  theory,  and  implementation  are  being  presented 
herein,  along  with  examples.  These  three  interface  elements  comprise  a  complete  set  of 
interftice  tools  for  global/local  analysis. 

1.  Introduction 

The  problem  of  connecting  dissimilar  meshes  at  a  common  interface  is  a  major  one  in 
finite  element  analysis.  One  method  of  connecting  these  dissimilar  meshes  is  to  use 
interface  elements. 

1.1  Applications 

The  interface  technology  can  be  divided  into  two  distinct  types.  The  first,  boundary 
technology,  connects  large  numbers  of  elements  along  a  sin^e  geometric  curve  or 
surface,  as  specified  by  the  user.  Primary  applications  facilitate:  global/local  analysis, 
vidiere  a  patch  of  elements  may  be  removed  from  the  global  -•  idel  and  replaced  a 
denser  patch  for  a  local  detail,  without  having  to  transition  to  i surrounding  area;  and 
component  assembly,  where  meshes  built  by  different  engineering  organi^ions  may  be 
connected,  such  as  a  vwng  to  the  fuselt^  of  an  urplane.  The  second,  interior  technology, 
connects  a  few  elements  together  in  a  mesh  transition  region,  wthout  user  specification. 
Primary  ^plications  enable:  automeshers,  which  may  be  required  to  transition  between 
large  and  small  elements  between  mesh  regions;  and  h-refinement,  where  subdivided 
elements  may  be  adjacent  to  undivided  elements  without  a  transition  area. 


1.2  Previous  Methods 

Much  work  has  been  done  to  resolve  the  element  interface  problem,  with  most  of 
the  efforts  concentrating  on  moving  the  nodes  or  writing  multi-point  constraint  (MPC) 
equations  on  the  interfaces.  However,  both  of  these  methods  have  difficulties  that  prevent 
them  from  being  practical  for  die  general  problem. 

13  Current  Method 

The  need  and  applications  for  reliable  interface  technology  are  great.  NASA 
Langley  Research  Center  has  developed  a  method  for  analyzing  structures  composed  of 
two  or  more  independently  modeled  substructures,  based  on  a  hybrid  variational 
formulation  with  Lagrange  multipliers,  and  applied  it  to  global/local  demonstration 
problems  for  one-dimensional  and  two-dimensiond  interfaces.  NASA  has  also  developed 
die  technology  for  a  solid-to-shell  transition  element  for  use  with  composites,  and  has 
combined  it  vrith  the  one-dimensional  interfrce  element 

Under  terms  of  a  cooperative  agreement  between  MSC  and  NASA,  MSC  has 
extended  and  implemented  this  interfrce  technology  into  MSC/NASTRAN  for  p-eiement 
edges  along  a  geometric  curve,  for  p-element  faces  over  a  geometric  surface,  and  for  solid 
p-element  faces  and  shell  p-element  edges.  This  agreement  is  part  of  NASA’s  continuing 
effort  to  transfer  technology  into  the  mainstream  of  industry  as  an  aid  in  developing 
competitiveness  in  the  worldwide  market. 
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Algorithms  that  find  good  partitionings  of  unstructured  and  irregular  graphs  are 
critical  for  the  efiicient  execution  of  a  wide  range  of  scientific  simulations  on  high 
performance  parallel  computers.  However,  for  many  types  of  computations  that 
are  now  commonplace,  such  as  adaptive  mesh,  multi-phase,  and  multi-ph3^ics  sim¬ 
ulations,  the  traditional  graph  partitioning  formulation  is  not  adequate.  For  these, 
more  generalized  graph  partitioning  formulations  and  algorithms  are  needed.  In  this 
talk,  we  describe  some  important  classes  of  scientific  simulations  that  require  new 
formulations  of  the  graph  partitioning  problem,  we  discuss  these  requirements,  and 
we  describe  new,  generalized  formulations  of  the  graph  partitioning  problem  as  well 
as  algorithms  for  solving  these  problem. 

For  large-scale  scientific  simulations,  the  computational  requirements  of  techniques 
relsdng  on  globally  refined  meshes  become  very  high,  especially  as  the  complexity  and 
size  of  the  problems  increase.  By  locally  refining  and  de-refining  the  mesh  to  capture 
flow-field  phenomena  of  interest,  adaptive  methods  make  standard  computational 
methods  more  cost  effective.  Similar  issues  also  exist  for  problems  in  which  the 
amount  of  computation  associated  with  each  mesh  element  changes  over  time.  For 
example,  in  particles-in-cells  methods  that  advect  particles  through  a  grid,  large 
temporal  and  spatial  variations  in  particle  density  can  introduce  substantial  lo^ 
imbalance.  In  both  of  these  types  of  applications,  it  is  necessary  to  dynamically 
load  balance  the  computations  as  the  simulation  progresses.  This  dynamic  load 
balancing  can  be  achieved  by  using  a  graph  partitioning  algorithm.  In  the  case  of 
adaptive  finite  element  methods,  the  graph  corresponds  to  the  mesh  obtained  after 
adaptation,  whereas  in  the  case  of  particles-in-cells,  the  graph  corresponds  to  the 

^This  work  was  supported  by  DOE  contract  number  LLNL  B347881,  by  Army  Research  Of¬ 
fice  contracts  DA/DAAG55-98- 1-0441  and  DA/DAAH04-95-1-0244,  by  Army  High  Performance 
Computing  Research  Center  cooperative  agreement  number  DAAH04-95-2-0003/contract  number 
DAAH04-95-C-0008,  the  content  of  which  does  not  necessarily  refiect  the  position  or  the  policy  of 
the  government,  and  no  official  endorsement  should  be  inferred.  Additional  support  was  provided 
by  the  IBM  Partnership  Award,  and  by  the  IBM  SUR  equipment  grant.  Access  to  computing 
facilities  was  provided  by  AHPCRC,  Minnesota  Supercomputer  Institute. 
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original  grid  with  properly  adjusted  vertex  weights  to  reflect  the  particle  density.  We 
will  refer  to  this  problem  as  adaptive  graph  partitioning  to  differentiate  it  from  the 
static  graph  partitioning  problem  that  arises  when  the  computations  remain  fixed. 
Adaptive  graph  partitioning  shares  most  of  the  requirements  and  characteristics  of 
static  graph  partitioning  but  also  adds  an  additional  objective.  That  is,  the  amount 
of  data  that  needs  to  be  redistributed  among  the  processors  in  order  to  balance  the 
computation  should  be  minimized. 

TVaditional  static  and  adaptive  graph  partitioning  formulations  are  not  general 
enough  to  ensure  the  elBBcient  execution  of  multi-phase  computations  on  high  per¬ 
formance  parallel  computers.  The  reason  is  that  these  applications  require  the  com¬ 
putation  of  partitionings  that  satisfy  more  than  one  balance  constraint,  while  tradi¬ 
tional  graph  partitioning  techniques  can  balance  a  single  constraint  only.  We  refer 
to  this  as  the  mtdti-constraint  graph  partitioning  pmblenu-As  in  the  case  of  a  single 
balance  constraint,  multi-constraint  partitioning  problems  can  be  either  static  or 
adaptive  in  nature.  Recently,  a  number  of  multi-constraint  graph  partitioners  have 
emerged  that  are  able  to  compute  partitionings  that  satisfy  a  number  of  balance 
constraints  while  also  minimizing  the  edge-cuts.  These  have  seen  a  great  deal  of 
success  in  meeting  the  partitioning  requirements  for  a  wide  range  of  multi-phase 
and  multi-physics  simulations.  However,  many  of  these  schemes  are  limited  for  two 
reasons,  (i)  They  are  serial  in  nature,  (ii)  They  address  only  the  static  multi¬ 
constraint  partitioning  problem  (and  not  the  adaptive  multi-constraint  partitioning 
problem). 

In  this  talk,  we  describe  formulations  of  the  multi-constraint  graph  partitioning 
problem  that  can  be  used  to  model  the  partitioning  requirements  of  a  wide  range 
of  scientific  simulations.  We  also  present  the  parallel  formulations  of  static  and 
adaptive  multi-constraint  graph  partitioning  algorithms  and  give  experimental  re¬ 
sults  from  128  processors  of  a  Cray  T3E.  We  show  that  our  parallel  algorithms  are 
able  to  efficiently  compute  partitionings  of  similar  quality  to  serial  multi-constraint 
algorithins,  can  scale  to  large  graphs,  and  are  very  fast.  For  example,  our  static 
multi-coristraint  partitioner  is  able  to  compute  a  three-constraint  128-way  parti¬ 
tioning  of  a  7.5  TtiillioTi  node  graph  in  about  7  seconds  on  128  processors  of  a  Cray 
T3E.  This  parallel  partitioner  is  used  as  the  key  component  of  an  adaptive  multi¬ 
constraint  algorithm.  We  show  that  the  adaptive  scheme  is  able  to  quickly  and 
robustly  balance  multi-constraint  partitionings,  while  also  minimizing  the  the  edge- 
cut.  Fi^hermore,  this  scheme  results  in  significantly  lower  data  redistribution  costs 
than  the  static  algorithm. 
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We  introduce  and  analyze  the  hp-version  of  the  Discontinuous  Galerkin  (DG)  time¬ 
stepping  method  for  non-linear  initial  value  ODEs  and  for  linear  parabolic  PDEs. 
New  a-priori  error  estimates  are  presented  which  are  completely  explicit  in  the  tirne 
steps,  in  the  approximation  orders,  and  in  the  regularity  of  the  exact  solution.  While 
these  estimates  allow  us  to  recover  the  optimal  convergence  rates  in  the  time  step  At, 
they  also  show  that  the  DG  method  converges  if  the  order  r  ->  oo  and  the  time  step 
is  kept  fixed.  We  are  able  to  prove  that  this  p-version  DG  approach  gives  spectral 
accuracy  for  solutions  with  smooth  time  dependence,  i.e.,  the  convergence  rates  are 
of  arbitrarily  high  algebraic  order.  Moreover,  for  entirely  analytic  solutions,  this 
p-version  convergence  is  of  exponential  order. 

Typically,  solutions  of  initial  value  problems  become  very  smooth  in  time  (in  many 
cases  even  analytic)  after  a  possibly  non-smooth  initial  phase  induced,  e.g.,  by  in¬ 
compatible  initial  data  or  by  piecewise  analytic  forcing  terms.  The  resolution  of 
such  solution  behavior  in  time  requires  appropriate  h-  and  p-refinement  towards 
the  singularity.  In  conjunction  with  geometric  time  steps  and  linearly  increasing 
approximation  orders,  we  show  the  hp-version  of  the  DG  method  can  approximate 
temporal  singularities  at  exponential  rates  of  convergence. 

A  complete  hp  discretization  in  time  and  space  is  discussed  for  linear  parabolic 
problems.  At  each  time  step  a  system  of  possibly  singularly  perturbed  reaction- 
dififiision  equations  is  solved.  If  the  hp-Finite  Element  Methods  for  these  spatial 
problems  take  into  account  certain  mesh  design  principles  using  anisotropic  and 
geometric  mesh  refinement  techniques,  exponential  rates  of  convergence  in  time  and 
space  can  be  achieved. 

Ail  theoretical  results  are  confirmed  and  illustrated  in  a  series  of  numerical  examples. 
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The  aim  of  our  contribution  is  an  efficient  method  which  is  also  parallelizable  pro¬ 
viding  the  following  features: 

1.  The  computation  of  boundary  displacements  and  boundary  tractions  with  high 
accuracy  either  on  the  whole  boundary  curve  (n  =  2)  or  boundary  surface 
(n  =  3)  or  on  an  a-priori  chosen  subregion  of  the  boundary  and 

2.  the  computation  of  displacements,  strains  and  stresses  near  and  up  to  the 
boundary  with  high  accuracy. 


Our  method  is  based  on  the  decomposition  idea  applied  to  rather  different  parts  of 
the  boundary  integral  equation  method: 

a)  the  decomposition  of  the  boundary  curve  or  boundary  surface  into  an  overlapping 
geometric  partition, 

b)  the  decomposition  of  the  trial  space  for  the  desired  quantities  into  a  regular  coarse 
grid  space  and  local  spaces  on  fine  grids;  and 

c)  the  additive  decomposition  of  the  boundary  integral  operators  involved  into  a 
standard  principal  part  and  a  smoothing  remainder. 

As  an  example  of  our  approach  we  consider  the  Lame  equations  of  linearized  isotropic 
homogeneous  elasticity, 

/xAu-l-(A-f-/i)graddivu  =  0  in  ft. 


where  the  Cauchy  data  on  the  boundary  P  =  dft  are  the  boundary  displacements 
and  the  boundary  tractions,  respectively,  i.e. 

ii|^  =  (f  and  A(div«)n  -f  fi  -H  J^n^gradu^j  =  if. 

Here  the  complete  stress  <r  on  the  boundary  is  for  two-dimensional  problems  given 
by 

ern|j,  =  ^  and  =  {if 'i^fi+ + 'tj  +  X{if 


where  near  to  the  boundary  we  use  Hadamard’s  natural  local  coordinates  in  nor¬ 
mal  and  tangential  directions  of  the  boundary,  respectivdy;  in  2-d  with  arclength 
parametrization  s  of  F  and  fi,  the  exterior  unit  normal,  t  the  unit  tangent  vector. 
For  3-d  problems,  on  the  boundary  surface  F ,  the  approach  is  principally  the  same. 

To  exemplify  our  method,  let  us  consider,  for  simplicity,  the  Dirichlet  problem  in 
plane  elasticity  where  we  have  to  solve  a  boundary  integral  equation  for  the  Cauchy 
datum  ipf  e.g.  the  integral  equation  of  the  first  kind, 

(5)  =F:=  (5/  +  K)0  on  F . 

In  a  first  step,  this  equation  is  solved  on  the  whole  boundary  where  we  use  coarse 
grid  trial  functions  on  a  regular  grid  associated  with  the  mesh  width  H.  In  order 
to  obtain  high  order  convergence,  we  use  the  method  proposed  in  [1],  [2],  where 
we  approximate  simultaneously  tangential  derivatives  of  V*  on  the  same  grid  which 
are  used  for  higher  order  recovery.  Since  the  coarse  grid  involves  only  a  relatively 
small  number  of  degrees  of  freedom,  the  coarse  grid  solution  approximating  low 
frequencies,  can  be  obtained  with  relatively  small  cost.  Next  we  use  the  geometric 
decomposition  together  with  an  associated  partition  of  the  unity  in  order  to  localize 
the  boundary  integral  equation.  Now  the  fine  grid  solution  on  a  window,  i.e.,  on 
one  of  the  subdomains  Fo  C  F,  can  be  solved  by  using  a  local  Galerkin  method.  In 
addition,  we  split  the  operator  into  a  simplified  principal  part  Vo  and  a  remainder 
that  has  a  kemeWhich  is  less  singular.  Then  we  solve  the  local  boundary  integral 
equation  for  u  €  H~^To)  on  the  fine  grid  where  the  simplified  operator,  i.e. 

(6)  {Vou,7^=={VoU}Uh,v)-{u{Vuh-F),^  for  all  v€lf"5(Fo) 

defines  the  corresponding  local  infiuence  matrix  whereas  the  modified  right-hand 
side  contains  all  the  global  information  in  terms  of  the  coarse  grid  solution,  the 
corresponding  residual  and  the  remainder  operator.  Again  we  use  simultaneous 
approximation  of  tangential  derivatives  and  recovery  for  improving  the  efficiency 
of  the  local  method.  For  smooth  boundaries  we  present  corresponding  asymptotic 
error  estimates  which  are  related  to  the  two-grid  finite  element  approximations  in 
[4].  For  nonsmooth  domains,  the  method  can  be  modified  appropriately  and  also 
performs  astonishingly  efficient  [3]. 
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Non-conformity  in  the  hp  version  can  involve  incompatibility  in  both  the  degrees 
and  the  meshes  between  adjoining  subdomams.-  In-this  regard,  the  mortar  finite 
element  techniques  have  become  very  popular  in  joining  together  such  incompati¬ 
ble  hp  discretizations.  In  this  talk,  we  introduce  some  variants  of  the  traditional 
mortar  finite  element  and  discuss  the  stability  of  these  methods  by  introducing  a 
new  measure,  that  takes  the  place  of  the  usual  inf-sup  stability  constant.  Our  nu¬ 
merical  results  show  optimality  of  the  resulting  non-conforming  method  for  various 
h,p  and  hp  discretizations,  including  the  case  of  exponential  hp  convergence  over 
geometric  meshes.  We  also  present  results  for  the  performance  of  these  methods 
when  the  stiffness  of  the  differential  operator  varies  in  different  subdomains.  Three- 
dimensional  considerations  suggest  that  the  variants  introduced  will  be  much  easier 
to  generalize  to  arbitrary  meshes  than  the  traditional  mortar  finite  element  method. 
Some  computational  results  for  the  stokes  problem  and  the  mixed  elasticity  problem 
discretized  by  the  hp  mortar  finite  element  method  will  also  be  presented. 
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^ile  the  Legendre-  or  Chebyshev-spectral  approximations  for  PDEs  in  bounded 
doTnains  have  achieved  great  success  and  popularity  in  recent  years  (see  e.g.  [3,  2]), 
spectral  approximations  for  PDEs  in  unbounded  domains  have  only  received  limited 
attention.  Pioneer  work  on  Laguerre  approximation  was  developed  in  Gottlieb  and 
Orszag  [3]  and  Maday,  Pemaud-Thomas  and  Vandeven  [5].  Most  of  the  previous 
numeric^  results  used  Laguerre  polynomials  which,  from  a  theoretical  point  of  view 
only  provide  meaningful  results  inside  small  intervals  since  the  error  estimate  is 
obtained  in  weighted  Sobolev  spaces  with  an  exponentially  decaying  weight,  and 
from  a  practical  point  of  view  are  not  suitable  for  practical  computations  due  to  the 
extremely  ill-conditioned  behaviors  of  the  Laguerre  polynomials  and  of  the  Laguerre- 
Gauss-Radau  quadrature  formula.  Furthermore,  it  was  pointed  out  by  Gottlieb  and 
Orszag  [3]  that  Laguerre  polynomials/functions  have  very  poor  resolution  properties 
when  compared  with  other  type  of  orthogonal  polynomials.  These  are  the  three 
main  factors  which  have  limited  the  interests  and  developments  on  using  Laguerre 
polynomials  for  problems  in  unbounded  domains. 

The  aim  of  this  work  is  to  develop  effective  remedies  for  the  aforementioned  diffi¬ 
culties  associated  with  Laguerre  approximations.  More  precisely,  we  will  propose 
uniformly  convergent  Galerkin  approximations  using  Laguerre  functions  for  elliptic 
equations  in  regular  unbounded  domains,  and  construct  stable  and  efficient  numer¬ 
ical  algorithms  for  their  practical  implementations.  Our  theoretical  error  estimates 
and  numerical  experiments  indicate  that 

•  from  both  the  theoretical  and  practical  points  of  view,  the  Laguerre  polyno¬ 
mials  are  generally  not  suitable  for  practical  implementations; 

•  with  a  proper  scaling,  the  Galerkin  method  using  Laguerre  functions  provides  a 
very  efficient  yet  very  simple  way  for  solving  problems  in  unbounded  domains; 

•  the  Laguerre-Galerkin  approximation  may  converge  exponentially  for  solutions 
which  delay  algebraicly  and  monotonically  at  infinity.  However,  only  algebraic 
convergence  is  possible  if  the  solution  decays  algebraicly  and  oscillates  at  in¬ 
finity. 
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In  summary,  as  demonstrated  in  this  paper  and  in  [4]  where  a  similar  Laguerre- 
Galerkin  method  is  successfully  used  for  the  approximation  of  some  nonlinear  partial 
differential  equations  in  semi-infinite  intervals,  we  believe  that  properly  formulated 
spectral  methods  using  Laguerre  functions  are  very  valuable  tools  and  viable  alter¬ 
natives  to  rational  polynomials  for  problems  in  unbounded  domains. 
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High-order  (p-version)  are  capable  of  exponential  rates  of  convergence  [1].  This 
increased  approximation  power  allows-the-use  ofcottse-meshes  in  udiich  relatively  few 
elements  cover  large  portions  of  the  computational  domain.  Preserving  the  convergence 
rate  for  problems  over  curved  domains  requires  that  a  curvilinear  mesh  geometry 
representation  be  used.  This  paper  is  concerned  with  the  task  of  generating  valid 
curvilinear  meshes  for  general  non-rtumifold  solid  models  as  defined  in  commercial 
geometric  modeling  systems. 


CAD  systems  are  becoming  the  mainstay  of  engineering  design  environment  The 
effective  application  of  engineering  analyses  requires  tiie  ability  to  directly  interact  udth 
these  models.  Over  the  past  several  years  the  technologies  have  been  developed  to 
siqjport  the  automation  of  low-order  finite  element  methods  for  objects  defined  in  a  solid 
modeling  system.  The  most  powerful  of  these  capabilities  directly  access  the  boundary 
representation  of  the  object  to  be  analyzed  and  employs  the  functionality  of  the  solid 
modeler  to  support  the  geometric  operations  needed  by  the  mesh  generator.  For  example 
the  MEGA  system  [2]  meshes  general  non-manifold  geometric  domains  through  a  set  of 
operators  that  have  been  directly  linked  with  the  geometry  engines  of  the  ACIS,  Parasolid 
and  Shapes  modeling  kernels  used  by  many  of  the  commercial  CAD  tystems.  From  the 
standpoint  of  siq)porting  engineering  aruilysis,  the  importance  of  non-manifold  solid 
models  is  tiiat  they  can  represent  objects  that  are  combinations  of  solids,  surfaces  and 
curves. 


This  paper  addresses  the  additional  capabilities  needed  to  generate  coarse  curved 
meshes  appropriate  for  p-version  tinite  element  methods.  The  methods  to  be  presented 
wll  rely  on  a  set  of  mesh  modification  operations  to  curve  those  mesh  entities  on  the 
boundaty  of  the  model  [3].  As  with  any  operations  that  modify  tiie  geometry  of  the  mesh, 
the  process  of  curving  mesh  entities  can  convert  valid  straight-sided  elements  to  invalid 
curved  elements.  Since  the  desired  meshes  can  be  coarse  with  respect  to  the  local  model 
curvature,  the  geometric  changes  associated  with  curving  the  boimdary  entities  can  be 
large.  Therefore,  the  algorithm  reqmred  to  successfully  create  the  curved  meshes  must 
include  a  complete  set  of  mesh  modification  operators  Aat  must  be  applied  in  a  carefully 
controlled  manner.  The  overall  algorithm  to  perform  these  modifications  is: 
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Preprocess  list  of  entities  to  be  modified  and  determine  target  geometry 

while  list  of  mesh  entities  not  empty,  loop  over  the  mesh  entities  in  the  list 

determine  allowed  move  that  maintains  a  limit  on  mesh  deformation  in  a 
single  step 

if  the  current  entity  can  accept  the  allowed  move 
move  to  that  location 

if  the  move  completes  the  move  to  the  target  geometry  remove  entity  from 
list 

else 

determine  the  first  mesh  entity  preventing  motion 
‘‘analyze”  current  situation  to  determine  “besT  local  modification 
perform  selected  local  modification 
determine  safe  move  distance  and  move 
determine  new  for  the  next  step 
end  if 
end  while 

The  operators  used  in  this  process  are  local  mesh  modifications  fruit  focus  on 
altering  the  mesh  entities  that  constrain  the  desired  mesh  motion  from  being  performed. 
The  basic  single  step  mesh  modifications  include  mesh  entity  swap,  collapse,  split  and 
geometry  change  operations.  Often,  it  is  not  possible  to  attain  the  desired  result  without 
chaining  a  small  number  of  these  operators  into  compound  operations. 


The  paper  mil  present  the  technical  details  of  a  geometry-based  mesh 
modification  procedure  for  the  generation  of  p-version  meshes.  Areas  to  be  addressed 
include  the  geometric  representation  of  the  mesh  entities,  the  mesh  modification  operators 
for  curved  meshes,  procedures  for  evaluating  the  geometric  properties  of  mesh  entities, 
and  the  details  of  the  mesh  curving  algorithm. 
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Trellis  is  a  geometry-based  analysis  fiameworic  for  finite  element  and  other  types  of 

numerical  analysis,  i^ch  is  based  on 

•  A  set  of  geometry-based  structures  which  can  support  direct  linkage  with  company 
CAD  information,  and  adaptivity  without  introducing  ^metric  ^proximation 
errors. 

•  A  careful  decomposition  of  die  geometry,  physics,  mathematical  model, 
discretization  and  numerical  methods  into  interacting  classes  that  siqiport 
multiphysics  coupling. 

•  Adaptive  control  of  the  numerical  solution  process. 

•  Simulation  automation  supported  by  the  MEGA  geometry-based  mesh  generation 
procedures. 


As  the  simulations  become  more  detailed,  the  computational  requirements  can  grow  too 
large  for  effective  solution  in  serial.  In  this  case  it  becomes  desirable  to  employ  parallel 
computing  techniques.  Ttellis  makes  use  of  the  Rensselaer  Partition  Model  (RPM)  as  tiie 
primary  tool  to  d^  with  the  parallel  data  management  issues  arising  fix)m  the  evolving 
nature  of  ad^tive  discretizations.  The  solution  of  the  linear  system  of  equations  makes 
use  of  genetic  programming  techniques  allowing  to  template  the  parallel  solution 
procedure  over  the  serial  implementatiotL 


This  prqrer  will  present  progress  to  date  on  the  development  of  the  geometry-based 
parallel  analysis  framework  Trellis.  Conceptually  Trellis  is  built  on  the  view  of  an 
analysis  as  a  transformation  between  three  levels  of  description.  The  highest  level 
description  is  that  of  the  physical  problem  which  is  posed  in  terms  of  physical  objects 
interacting  with  their  environment  Since  the  goal  of  the  analysis  is  to  obtain  reliable 
estimates  of  the  response  of  the  system  the  second  level  is  a  mathematical  problem 
description  that  introduces  some  level  of  idealization,  which  also  needs  to  be  controlled 
to  yield  the  desired  accuracy.  The  third  level  is  the  numerical  discretization  constructed 
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from  a  mathematical  problem  that  involves  another  set  of  idealizations  that  also  need  to 
be  controlled. 


Trellis  represents  a  careful  implementation  of  the  overall  structure  of  a  numerical  analysis 
process  into  a  framewoik.  It  starts  at  die  mathematical  problem  description  providing 
classes  diat  predefine  the  design  so  that  the  application  programmer  can  concentrate  on 
the  specifics  of  the  application.  Key  abstractions  in  Trellis  are  the  analysis  class  that  is 
being  used  to  transform  a  mathematical  problem  into  a  discrete  form.  The  discrete  form  is 
represented  by  the  discrete  system  class  that  contains  the  overall  system  represented  as 
contributions  from  system  contributors  (stiffness  contributors,  force  contributors,  and 
constraints).  At  diis  point  the  problem  is  described  in  an  analysis-independent  maimer. 
Constructing  the  discrete  system  in  parallel  is  based  on  utilizing  RPM,  which  provides  a 
notion  of  partitions,  and  partition  toundaries  describing  the  distributed  entities  of  the 
discretization. 


The  solution  of  the  linear  systems  of  equations  in  parallel  is  based  on  generic 
programming  techniques.  The  solver  components  are  written  to  conform  to  specified 
requirements  so  that  the  actual  parallel  solver  can  be  templated  over  die  serial 
implementation  of  matrices,  vectors,  and  operations  on  them,  as  well  as  a  connectivity 
class  describing  the  representation  of  the  p^lel  distribution.  This  allows  die  extension 
of  the  solver  part  wito  new  matrix  and/or  vector  storage  formats,  without  the  need  to 
specify  parallel  implementations  for  the  algorithms. 
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The  development  of  unstructured  spectral/hp  elements  using  elemental  domains 
such  as  tetrahedrons  and  triangles  has  encouraged  the  use  of  automated  mesh  gen¬ 
eration.  However  most  mesh  generators  have  been  developed  for  use  with  linear 
elements. 

Although  all  the  problems  associated  with  mesh  generation  for  linear  elements  are 
inherent  in  mesh  generation  for  high  order  elements  the  converse  is  not  true.  For 
example,  the  introduction  of  surface  curvature  in  the  elemental  representation  can 
lead  to  self  intersecting  elements  which  will  have  a  vanishing  or  negative  Jacobian 
thus  a  singular  isoparametric  mapping.  This  problem  is  also  exaggerated  by  the  fact 
that  when  using  high  order  elements  we  typically  want  a  much  coarser  discretisation. 

Even  when  the  surface  is  not  curved  the  mapping  from  the  parametric  surface  repre¬ 
sentation  to  the  physical  coordinates  can  generate  problems  as  illustrated  in  figure  2. 
In  this  example  we  have  represented  the  planar  surfaces  of  a  cube  by  an  equispaced 
mesh  of  discrete  points  and  a  quadratic  distribution  of  discrete  points.  In  the  case 
of  the  equispaced  definition  the  mapping  to  the  parametric  space  of  the  representa¬ 
tion  of  the  surface  is  isometric  and  we  obtain  the  surface  distribution  of  elemental 
points  shown  on  the  left.  However  in  the  case  where  we  use  a  surface  mapping  with 
quadratic  distribution  of  points  we  obtain  the  distorted  surface  reconstruction  of 
the  cube  boundary  shown  on  the  right. 

These  problems  have  been  highlighted  from  investigations  in  the  field  of  compu¬ 
tational  haemodynamics  which  involve  the  generation  of  computational  meshes  of 
arterial  bypass  grafts.  The  paper  will  address  the  problems  associated  with  sur¬ 
face  modelling  of  highly  complex  geometries  and  discuss  strategies  to  make  the  hi^ 
order  mesh  generation  and  adaption  more  accurate  and  robust. 
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(a)  (b) 

Figure  2:  Deformation  of  surface  element  due  to  the  mapping  from  the  parametric 
space  of  the  surface  representation.  The  mapping  in  (a)  is  isometric  and  the  mapping 
in  (b)  is  not. 
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The  application  of  three-dimensional  computational  modelling  to  study  the  human 


Figure  3:  Wall  shear  stress  distribution  over  a  reconstruction  of  a  coronary  bypass 

arterial  system  has  become  widespread  over  the  last  decade.  This  area  presents  a 
series  of  ^allenges  in  flow  modelling  due  to  pulsatile  nature  of  the  flow  which  is  typ¬ 
ically  at  a  moderate  Reynolds  number  regime  where  both  viscous  and  inertial  effects 
are  important.  Further  we  recognise  that  the  flow  geometries  are  very  complex,  the 
artery  walls  are  distensible  and  the  blood  behaviour  is  non-Nevirtonian. 

The  good  accuracy  and  phase  properties  of  high  order  algorithms  has  made  them 
particularly  suitable  for  scientific  studies  of  flow  physics  and  therefore  for  under¬ 
standing  haemodynamics.  Furthermore  there  is  a  strong  body  of  evidence  to  show 
that  disease  processes  such  as  myointimal  hyperplasia  and  arteriosclerosis  are  re¬ 
lated  to  wall  shear  stress.  Therefore  there  is  a  need  to  accurately  model  both  the 
primitive  variables,  v  &  p,  and  the  derivative  quantities  associated  with  wall  shear 
stress. 
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Figure  3  shows  a  computational  reconstruction  from  an  MRl  image  of  an  coronary 
bypass  graft  and  the  computed  wall  shear  stress  distribution.  The  paper  will  discuss 
how  high  order  unstructured  spectral//ip  element  methods  are  being  used  to  model 
flows  in  arterial  bypass  grafts  and  the  possible  implication  for  interventional  surgery. 
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We  consider  screen  problems  (with  a  hypersingular  or  weakly  sin^ar  integral  equa¬ 
tion  on  an  open  surface  T)  as  model  problems.  For  the  hypersingular  equation  the 
preconditioner  is  based  on  a  three-level  decomposition  of  the  underlying  ansatz 
space,  the  levels  being  piecewise  bilinear  functions  on  a  coarse  grid,  piecewise  bilin¬ 
ear  fimctions  on  a  fine  grid,  and  piecewise  polynomials  of  high  degree  on  the  fine  grid. 
We  prove  that  the  condition  number  of  the  preconditioned  linear  sjrstem  is  bounded 
by  maxj(l  -I-  log  where  Hj  is  the  diameter  of  an  element  of  the  coarse  grid, 
hj  is  the  size  of  the  elements  of  the  fine  grid  on  Tj  and  pj  is  the  maximum  of  the 
polynomial  degrees  used  in  F^-.  For  the  weakly  singular  integral  equation,  where  no 
continuity  of  test  and  trial  functions  across  the  element  boundaries  has  to  be  en¬ 
forced,  the  method  works  for  non-uniform  degree  distributions  as  well.  We  comment 
alsn  on  various  adaptive  algorithms.  Numerical  results  supporting  our  theory  are 
reported. 
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Experimental  determination  of  loads  that  can  cause  buckling  and  other  fsdlues  is 
both  expensive  and  unreliable.  Mathematical  models  used  in  current  engineering 
practice  are,  on  the  other  hand,  based  on  dimensionally  reduced  descriptions  of  an 
elastic  body.  The  assumptions  necessary  for  such  dimensional  reduction  to  hold  can 
often  result  in  a  large,  unknown  modeling  error  when  non-zero  initial  stress  sta.tes 
are  present,  leading  to  inaccurate  failure  predictions. 

A  linearized  model  for  buckling  and  stress-stiffening  by  Szabo  and  Kiralyfalvi  [1] 
has  been  recently  implemented  in  the  hp  code  STRESSCHECK.  This  model  does 
buckling  analysis  for  the  fully  three-dimensional  problem  at  hand,  rather  than  some 
asymptotic  (dimensionally  reduced)  limit.  It  finds  the  smallest  positive  multiple  A 
of  an  existing  (pre-budding)  stress  state  (Tq  that  will  result  in  buckling.  The  use 
of  the  hp  method  enables  solutions  over  singular  domains  to  be  well  approximated, 
and  ensures  that  no  loddng  takes  place  even  when  the  domain  is  very  thin. 

However,  a  potentially  serious  danger  of  the  method  is  that  it  characterizes  A  as  the 
lowest  positive  spectral  value  of  a  problem  of  the  form 

A  -  AB  =  0 

where  A  and  B  are  both  second-order  differential  operators.  The  finite  element 
approximation  of  such  (non-compact)  problems  could  (at  least  theoretically)  be 
very  ill-behaved,  due  to  the  presence  of  spurious  approximate  eigenvalues,  which 
can  completely  pollute  the  results. 

We  show  that  (1)  spurious  eigenvalues  are  absent  for  problems  of  engineering  interest 
for  which  the  pre-buckling  stress  cq  is  bounded  (2)  spurious  eigenvalues  are  always 
present  for  problems  where  oo  is  unbounded.  For  the  latter  case,  we  demonstrate 
how  the  reliability  of  the  computations  can  still  be  assessed,  using  the  eigenvectors 
(the  buckling  shapes).  We  describe  our  results  below. 

1.  Absence  of  spurious  eigenvalues  for  bounded  <ro 

By  considering  the  asymptotics  of  a  plate  of  thickness  d  ,  we  show  that  no  polluting 
eigenvalues  will  be  observed  in  an  interval  \C\d~^y  C72d“^]  about  the  origin.  Hence, 


spectral  pollution  will  not  occur  for  the  eigenvalues  of  interest,  provided  the  domain 
is  thin.  (The  thin  case  is  the  case  of  engineering  interest.) 

Again  considering  the  asymptotics  of  a  plate,  we  investigate  the  effect  of  various 
initial  stress  states.  We  show  that  the  only  case  where  spurious  eigenvalues  could 
occur  is  the  one  where  the  initial  stress  is  of  pure  bending  type.  This  case  is, 
however,  uninteresting  from  the  point  of  view  of  applications.  For  all  initial  stresses 
that  occur  in  applications  (these  can  be  decomposed  into  a  bending  and  a  non-zero 
membrane  stress),  spectral  pollution  does  not  occur  in  the  region  of  interest. 

Our  asymptotic  analysis  has  also  established  that  boundary  layers  occurring  in  the 
solution  do  not  lead  to  adverse  affects  in  the  approximation  of  the  desired  eigenval¬ 
ues. 

Finally,  we  characterize  the  initial  stresses  that  lead  to  pollution-free  eigenvalue 
approximations  when  parts  of  the  domain  are  not  thin. 

2.  Reliability  of  computations  for  unbounded  <to 

For  domains  with  comers,  or  domains  with  a  change  in  the  type  of  boundary  condi¬ 
tions,  it  is  well-known  that  the  stresses  can  have  r”®  singularities.  We  show  that  in 
this  case,  the  value  0  always  lies  in  the  continuous  spectrum  of  the  problem.  This 
causes  the  numerical  approximations  (which  appear  as  a  cluster  of  spurious  eigen¬ 
values)  to  converge  to  0.  In  applications,  however,  one  is  interested  in  the  lowest 
eigenvalue  of  a  family  related  to  the  spectrum  of  a  limiting  problem  (the  d  =  0  case). 
We  characterize  the  values  of  a  for  which  this  limiting  spectrum  is  still  well-defined, 
and  well-behaved. 

We  show  computationally  determined  buckling  shapes  associated  with  spurious  and 
non-spurious  modes.  It  is  seen  that  for  spurious  modes,  the  eigenfunction  is  almost 
constant  over  the  whole  domain,  except  for  the  region  surrounding  the  point  of 
singularity.  These  can  therefore  be  easily  distinguished  from  non-spurious  modes, 
where  the  eigenfunction  is  ‘global’  as  opposed  to  ‘local.’  Moreover,  although  the 
spurious  eigenvalues  converge  exponentially  to  0,  the  limiting  eigenvalues  of  interest 
converge  even  faster,  and  hence  can  be  recovered,  as  long  as  they  are  smaller  than 
the  spurious  eigenvalues. 

Two  important  facts  that  emerge  from  our  investigation  are:  (1)  The  degree  of 
localness  of  the  eigenfunctions  can  be  an  effective  tool  in  assessing  the  reliability  of 
the  computations  (2)  There  is  often  an  optimal  level  of  layer  refinement,  over  which 
spurious  eigenvalues  can  arise.  Our  results  therefore  indicate  that  subject  to  the 
above  limitations,  the  hp  method  in  [1]  is  robust  for  engineering  applications. 
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For  the  past  few  decades  more  and  more  industries  are  using  computer  simulations  to 
establish  the  response  of  the  structure  to  loads  associated  with  operational  conditions.  The 
advantage  of  the  computer  simulation  over  a  more  traditional  ^proach  (make  and  break) 
is  obvious:  companies  save  millions  of  dollars  >^en  costly  experiments  are  replaced  by 
analyses.  However,  computer  simulation  of  die  behavior  of  a  real  structure  is  not  a  trivid 
task.  The  finite  element  method,  provided  in  a  number  of  commercial  codes,  is  a 
commonly  used  numerical  simulation  tool.  Although  the  finite  element  tool  is  the  most 
popular  simulation  tool  used  in  industry,  the  capabilities  and  limitations  of  the  tool  are 
seldom  well  understood  by  the  tools*  users.  It  is  "NEWS"  to  many  users  that  the  finite 
element  method  is  just  a  numerical  approximation  to  the  solution  to  a  particular  set  of 
partial  differential  equations  with  prescribed  boundary  conditions — in  structural  analysis, 
of  course  these  are  the  equations  of  equilibrium  (motion)  for  a  single  body.  Real  world 
structures,  however,  consist  of  many  parts,  and  die  task  is  to  assess  the  behavior  of  the 
entire  assembly.  Mathematical  modeling  of  the  connection  between  different  parts 
becomes  an  essential  part  of  many  simulations.  A  variety  of  contact  algorithms  are 
implemented  in  the  commercial  codes  that  perform  this  task.  It  is  ho  -  ver  questionable 
whether  the  avmlable  tools  provide  reliable  results  which  can  b;  used  to  closely 
approximate  the  r^ponse  of  the  structure.  In  diis  paper  finite  element  models  of  a 
structural  assembly  diat  is  a  part  of  a  combustor  area  in  an  industrial  turbogenerator  are 
aiuilyzed.  Figure  1.  The  connectivity  between  different  parts  in  this  model  is  often 
simulated  with  contact  or  gap  elements.  In  this  piq)er,  a  small  subset  of  this  larger  more 
complicated  problem  will  be  analyzed.  Results  obtained  using  /i-version  commercial  code 
ANSYS  are  compared  with  results  using  the  /7-version  code  StressCheck.  The  goals  of 
computation  and  die  quality  of  the  results  are  discussed. 
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Figure  1  Combustor.  Top-Solid  Model.  Bottom-Mesh 
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Introduction 

We  propose  to  demonstrate  the  application  of  the  discontinuous  Galerkin  method 
(DGM)  to  Maxwell’s  equations  using  new  developments  for  nodal  unstructured  spec¬ 
tral  elements  in  two  and  three  dimensions.  Specifically  we  will  discretize  the  evo¬ 
lution  equations  for  the  magnetic  field  and  electric  field  vectors  in  a  medium  of, 
possibly  sharply,  variable  permeability  and  permittivity.  This  is  the  typical  physi¬ 
cal  description  for  the  electromagnetic  fields  emitted  by  a  mobile  phone  and  their 
transmission  in  the  air  and  user’s  skull  cavity.  There  are  currently  significant  de¬ 
bates  about  the  possible  harm  caused  by  prolonged  use  of  mobile  phones.  We  wish 
to  demonstrate  that  using  DGM  on  unstructured  grids  will  be  able  to  give  accu¬ 
rate  results  for  high  frequency  electromagnetic  waves  passing  through  such  sharply 
heterogeneous  materials  in  complex  geometries. 

Proposed  Numerical  Studies 

The  examples  we  will  use  are  scattering  by  perfectly  conducting  bodies,  scattering 
by  dielectrics,  and  scattering  by  sharp  points  (i.e.  a  plane  wave  incident  on  a  cone). 
These  last  two  are  most  diallenging  for  spectral  methods,  where  the  non-smoothness 
in  the  equations  and  boimdary  conditions  can  cause  0(1)  errors.  We  will  show  that 
the  DGM  handles  these  cases  eflaciently  and  robustly.  We  will  also  show  preliminary 
results  for  scattering  from  an  F15  geometry. 

Example:  There  is  an  exact  Mie  solution  for  scattering  of  incident  plane  waves  b; 
a  perfectly  conducting,  infinitely  long  cylinder.  The  exact  solution  for  the  electric 
field  component  is: 

EripA)  = -E,  t 

n=-oo  (Ka) 

in  terms  of  (p,  <f>)  polar  coordinates. 

We  have  already  run  this  simple  simulation  in  a  domain  of  approx  16  diameters, 
with  a  radial  ABC.  We  show  in  figure  4  the  numerical  solution  at  t  =  30  using  the 
DGM  approach.  This  is  visually  indistinguishable  from  the  exact  solution.  However 
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in  the  test  of  L2  and  Loo  error  (not  including  the  ABC  region)  we  see  that  the  error 
saturates  at  a  level  of  10“^.  This  saturation  is  due  to  reflections  from  the  ABC  and 
in  the  proposed  paper  we  intend  to  investigate  how  to  modify  this  and  alternative 
ABCs  to  reduce  back  scattering  from  far-field  boundaries.  We  will  demonstrate 
alternative  results  using  a  PML  (perfectly  matched  layer)  formulation  in  order  to 
improve  these  results. 


Figure  4;  Left:  DGM  solution  with  ABC  and  p=10  at  t=30  for  scattering  of  a  plane 
wave  incident  on  a  perfectly  conducting  cylinder.  Right:  Convergence  plot  for  L2 
and  Loo  error  for  the  DGM  approximation  compared  with  the  Mie  solution  at  t  =  30 
as  a  function  of  polynomial  order. 
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The  classical  formulation  of  the  spectral- elment- method  based  on  the  technique 
of  Maday  and  Patera  depends  on  the  choice  of  fimctional  space,  grid  points,  and 
method  of  integration.  On  the  quadrilateral,  the  preferred  space  is  a  tensor  product 
basis  of  Legendre  polynomials, 

(7)  Va  =  {span  (x"y’")|0  <m<N;0<n<  N}. 

The  grid  points  are  chosen  as  the  Gauss-Lobatto-Legendre  points,  which  lead  to 
an  accurate,  well-behaved  interpolation,  and  subsequently  well-behaved  derivatives. 
The  choice  of  points  is  also  important  for  coupling  elements  together  easily.  By 
choosing  a  quadrature  slightly  less  accurate  than  the  optimal  Gaussian,  one  achieves 
a  computationally  eflScient  diagonal  global  mass  matrix  for  explicit  methods. 

On  the  triangle  the  preferred  space  is, 

(8)  V\  =  {span(a:"y”*)|0  <m,n<K\m  +  n<  K). 

The  space  was  also  used  by  Dubiner  who  reintroduced  an  orthogonal  basis  on 
the  triangle,  which  was  first  presented  by  Koomwinder.  Dubiner  also  introduced  a 
warped  tensor  product  quadrature  set  for  the  triangle.  However,  this  quadrature 
set  is  oversampled  which  leads  to  a  dense  global  mass  matrix.  This  is  much  more 
expensive  than  for  the  cleverly  implemented  explicit  method  on  quadrilaterals. 

Triangle  Functional  Space 


On  the  quadrilateral  the  favorite  basis  functions  are  eigenfunctions  of  a  singular 
Sturm-Liouville  problem,  therefore,  by  the  standard  integration  by  parts  arguments, 
the  convergence  rate  for  the  approximation  of  smooth  functions  depends  on  the 
regularity  of  the  function  being  approximated  and  not  on  an  special  conditions  at 
the  element  boundaris.  It  has  been  known  for  quite  some  time  that  similar  problems 
exist  for  the  triangle.  We  now  know  the  following  about  simplices: 
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1.  We  now  know  the  Koornwinder  polynomials  are  the  eigenfunctions  of  a  singu¬ 
lar  Stunn-Liouville  problem  on  the  simplex.  Therefore,  the  convergence  rate 
for  the  approximation  of  smooth  functions  depends  on  the  regularity  of  the 
function  being  approximated  and  not  on  an  special  conditions  at  the  element 
boundaries  [4]. 

2.  Unlike  the  case  on  the  line,  we  now  know  the  eigenvalues  come  in  eigenspaces. 
The  commonly  used  triangle  truncation  is  a  union  of  these  eigenspaces  for 
all  eigenvalues  less  than  a  constant.  This  tnmcation  is  invariant  under  the 
natural  symmetry  group  of  the  triangle  [4]. 

3.  We  now  know  Jackson  and  Berstien  type  inequalties  for  these  spaces  [1]. 

Triangle  Interpcdation- Points 

There  has  been  recent  work  on  computing  optimal  interpolation  points  on  the  sim¬ 
plex.  Each  set  of  points  has  interesting  properties,  but  we  prefer  Fekete  points,  as 
in  [2],  for  the  following  reasons: 

1.  On  the  [—1, 1]  interval,  Fekete  points  are  the  Gauss-Lobatto  points. 

2.  On  the  square,  Fekete  points  are  the  tensor  product  of  Gauss-Lobatto  points. 
Thus,  the  conventional  spectral  element  method  can  also  considered  to  be  a 
Fekete  point  method. 

3.  Under  suitable  assumptions  we  can  show  that  the  Fekete  points  along  each  edge 
of  the  triangle  are  Gauss-Lobatto  points  making  the  triangles  and  quadrilat¬ 
erals  naturally  conform. 


Method 

We  introduced  the  Fekete  point  spectral  element  method  in  [3].  There  we  showed 
spectral  convergence  for  linear  advection  problems.  In  this  work  we  answer  some 
key  remaining  questions  about  the  method; 

1.  Will  the  points  integrate  twice  the  space  like  Gauss-Lobatto  points  on  the 
quadrilateral?  Does  the  integration  preserve  integration  by  parts? 

2.  Is  the  method  stable? 

3.  Is  the  method  as  efficient  as  that  on  the  quadrilaterals? 

4.  Does  the  method  perform  as  well  for  non-linear  problems? 
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The  solution  to  the  Laplace  operator  in  three-dimensional  domains  in  the  vicinity 
of  straight  edges  is  presented  as  an  asymptotic  expansion  involving  eigen-pairs  with 
their  coefficients  cailed  edge  flux  intensity  functions  (EFIFs).  The  eigen-pa^rs  are 
identical  to  their  two-dimensional  counterparts  over  a  plane  perpendicular  to  the 
edge.  The  general  form  of  the  asymptotic  expansion  is: 

u{r,0,X3)  =  535j(r,d,a:3) 

i=l 

=  ^  ai{xz)r**‘si{9)  +  Ciid^ai{xz)r^*'^\{e)  +  cadzOiMr***'^  +  •■■ 

t=t 

where  xz  is  the  coordinate  along  the  edge,  and  r,  6  are  cylindrical  coordinates  in  a 
plane  perpendicular  to  the  edge,  aj+i  >  are  called  edge  eigen-values,  ai(x3)  are 
analytic  in  0:3  and  are  denoted  by  edge  flux  intensity  functions  {EFIFs).  3i{0)  are 
analytic  in  0,  called  edge  eigen-function  and  Cij  are  given  constants. 

Extraction  of  EFIFs,  cannot  be  obtained  in  a  straightforward  manner  over  this 
two-dimensional  plane.  Consequently,  a  special  method  based  on  projection  and 
Richardson  extrapolation  is  presented  for  point-wise  extraction  of  EFIFs  &om  p- 
finite  element  solutions.  The  mathematical  analysis  is  demonstrated  by  numerical 
e^erimentation.  A  similar  but  more  efficient  method,  based  on  "energy  projection", 
for  extracting  EFIFs  is  proposed.  The  details  are  provided  in  [1]. 
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Viscoelastic  behavior  has  been  become  of  great  importance  with  modem  develop¬ 
ments  in  technology,  such  as  the  demands  on  hi^  temperature  behavior  of  materials 
for  gas  turbines,  jet  engines  and  power  plants,  and  also  with  extensive  studying  of 
biomaterials.  However,  the  comprehensive  understanding  of  the  mechanics  of  vis¬ 
coelastic  problems  is  still  hindered  by  a  lack  of  reliable  numerical  simulation  methods 
and  experimental  data  for  the  characterization  of  material  properties.  So  far,  most 
of  the  procedures  implemented  in  various  finite,  element  software  are  based  on  the 
traditional  h-version  of  the  finite  element  method.  These  procedures  usually  do  not 
provide  for  control  of  the  solution  quality,  making  it  difficult  to  estimate  and  con¬ 
trol  discretization  errors.  Proper  control  of  discretization  errors  is  cmcial  for  the 
selection  of  the  proper  constitutive  laws  and  the  determination  of  the  coefficients  in 
the  constitutive  laws. 

Our  goal  is  to  develop  a  stable  and  reliable  algorithmic  procedure  for  the  numerical 
i^ulation  of  viscoelastic  solids  which  allows  estimation  and  control  of  the  errors 
of  discretization,  therefore  the  material  properties  can  be  obtained  efficiently  iising 
optimization  algorithms  that  minimize  the  difference  between  experimental  measure¬ 
ments  and  finite  element  predictions.  Towards  this  end  we  have  developed  a  new 
algorithmic  procedure  to  solve  viscoelastic  problems  using  the  mixed  formulation. 
The  control  of  spatial  discretization  errors  is  based  on  the  hp-version  of  the  finite 
element  method  in  which  the  displacement  and  stress  fields  are  approximated  sep¬ 
arately.  The  control  of  time  integration  errors  is  based  on  finite  difference  method. 
The  algorithm  has  been  implemented  into  a  research  code  for  planar  and  axisym- 
metric  elements.  Both  linear  and  nonlinear  constitutive  laws  are  incorporated.  For 
linear  viscoelastic  materials,  the  general  differential  constitutive  equation  is  used. 
For  nonlinear  viscoelastic  materials,  several  empirical  creep  laws  are  used.  Small  dis¬ 
placement  and  small  strain  conditions  are  assumed.  Numerical  results  are  presented 
for  benchmark  problems  which  show  the  robustness  of  the  formulation. 

Numerical  testing  has  been  performed  with  the  objective  to  investigate  the  locking- 
free  properly  of  the  mixed  formulation  for  nearly  incompressible  materials.  The 
results  show  that  the  method  is  free  from  locking  and  performs  well  for  incompress¬ 
ible  materials. 
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